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I Introduction 



The Standard Model (SM) (Weinberg, 1967; Salam, 1969; Glashow et ai, 1970) of strong 
and electroweak interactions has been tested with precise experiments at the present coUid- 
ers (for reviews see, e.g., AltareUi, 1989; Marciano, 1991, 1993a,b; Bethke, 1992; Bethke 
and Pilchcr, 1992; R. K. Elhs, 1992; Langackcr, Luo and Mann, 1992; Brodsky, 1993; 
Sirhn, 1993a,b; Kniehl, 1994a; Soper, 1995). However, the decisive confirmation of the 
SM or its modification is still ahead. One awaits with great interest precise measurements 
from LEP, SLC, HERA, Fermilab Tevatron, etc. 

The progress in a very precise experiment requires adequate progress in the developing 
of calculational methods and performing the theoretical computations of various observ- 
ables. This can be consistently done within perturbation theory for processes with large 
momentum transfer. Nowadays, the standard way to evaluate experimentally measur- 
able quantities from first principles of the theory is to use perturbation methods. Lattice 
calculations provide an alternative method. 

The main goal of this paper is to review some of the recent achievements in methods 
of high order analytical perturbative calculations of a wide class of observable quanti- 
ties. These quantities are total cross sections, decay widths and structure functions in 
deep inelastic processes, several key theoretical quantities, such as renormalization group 
functions, renormalization constants, Wilson coefficient functions, etc. We will present a 
simplified description of some of the recent calculations. 

A decisive role in the construction of the SM has been played by experimental studies 
of so called inclusive processes, in particular deep inelastic lepton-hadron processes like 
e'^e~ — > hadrons, deep inelastic e, /i and v -scattering, etc. The discovery of scaling of deep 
inelastic structure functions (Bjorken, 1968, 1969; Yang, 1969) led to the parton model 
(Feynman, 1969, 1972; sec also DrcU and Yan, 1971). The explanation of the observed 
scaling properties has been given by Matveev, Muradyan and Tavkhclidze (1970, 1972), 
using the universal principle of automodelity and dimensional analysis. The quark count- 
ing formulae, allowing one to obtain the high energy asymptotic behavior for cross sections 
and hadron form factors at large momentum transfers, have been derived by Brodsky and 
Farrar (1973), and Matveev, Muradian and Tavkhelidze (1973). The discovery of asymp- 
totic freedom in nonabelian gauge field models (Gross and Wilczek, 1973; Politzer, 1973) 
together with the conception of spin half, fractionally electric charged fundamental con- 
stituents of hadrons - quarks (Gell-Mann, 1964; Zweig, 1964) with an additional quantum 
number color (Bogolyubov, Struminsky and Tavkhclidze, 1965; Tavkhelidze, 1965; Han 
and Nambu, 1965; also Miyamoto, 1965; Grccnberg, 1964), interacting via eight massless, 
non-abelian, spin 1, self interacting gauge fields - gluons, led to the creation of Quantum 
Chromodynamics (QCD) (Fritzsch, Gell-Mann and Leutwyler, 1973) - the present theory 
of strong interactions. For an introductory review of QCD see, e.g., Marciano and Pagcls 
(1978), and a short historical review has been given recently by Tavkhelidze (1994). QCD 
is based on a local SUc(3) symmetry group, which implies the minimal locally gauge in- 
variant Lagrangian density of the model. QCD is a renormalizable quantum field model. 
There exists well defined rules for removing of ultraviolet divergences from S-matrix am- 
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plitudes at each order of the interaction couphng constant. After the renormahzation, the 
calculated physical quantities are free of ultraviolet regularization parameters. The prob- 
lem of renormalizability of non-abelian gauge theories has been considered since the early 
60's (Feynman, 1963; deWitt, 1967; Mandelstam, 1968, etc.). After the Lorentz covariant 
quantization of gauge fields, based on the path integral approach (Faddeev and Popov, 
1967; for a textbook, see also Faddeev and Slavnov, 1980), the proof of renormalizability 
was given ('t Hooft, 1971; Lcc and Zinn- Justin, 1972, 1973). Besides the short distance 
effects, in QCD one has to deal with infrared divergences associated with long distance 
infinities. In other words, in addition to the large parameter (large momentum transferred 
- Q^), there are small parameters such as, for instance, hadron mass - m or momenta of 
some of the participating particles, and in the calculation one faces senseless large loga- 
rithmic contributions ~ \ogm? / . The infrared divergence problem was considered long 
ago for QED (Bloch and Nordsieck, 1937; Yennie, Frautschi and Suura, 1961). A modern 
treatment of this problem in the SM is based on the operator product expansion technique 
(Wilson, 1969) and factorization theorems. This, in some cases, in particular, for deep in- 
elastic processes, allows one to factorize the large and small distance contributions (Libby 
and Sterman, 1978; Mueller, 1978; Efremov and Radyushkin, 1980a,b; Radyushkin, 1983; 
CoUins and Soper, 1987; Collins, Soper and Sterman, 1983, 1984, 1985, 1989 and refer- 
ences therein). The concrete prescriptions of dealing with infrared divergent Feynman 
integrals have been given by Vladimirov (1978, 1980), Pivovarov and Tkachov (1988), 
Tkachov (1991, 1993), Chetyrkin and Smirnov (1985). For earlier references, see the work 
by Tkachov (1993). For a textbook, see CoUins (1984). The interference of long and short 
distance effects is still problematic in actual higher order calculations. 

The group character of renormalizing transformations in quantum field theory was 
first discovered by Stueckelberg and Peterman (1953), and Gell-Mann and Low (1954) 
have applied it to study the ultraviolet asymptotics of Green's functions in spinor elec- 
trodynamics. The mathematical formalism of the renormahzation group has been worked 
out by Bogolyubov and Shirkov (1955, 1956a,b), and Bogolyubov and Parasyuk (1955a,b, 
1956, 1957) have introduced the i?-operation for subtracting ultraviolet divergent contri- 
butions recursively on the level of loop Feynman diagrams. The renormahzation group 
and /^-operation techniques are the crucial tools in any perturbative calculation within 
the Standard Model (see the textbook by Bogolyubov and Shirkov, 1980). For a historical 
review on renormahzation group sec, for example, Shirkov (1992; also Peterman, 1979) 
and references therein. For a textbook on the modern renormahzation theory and the 
references, see Collins (1984). 

The property of asymptotic freedom, the method of renormahzation group and fac- 
torization theorems are the basis of the present perturbative QCD. In order to relate 
perturbative QCD to measurable quantities, it was necessary along with the renormahza- 
tion group, dispersion relations, operator product expansion techniques and factorization, 
to develop a technique for evaluation of loop Feynman diagrams. Indeed, in each order 
of perturbation theory, contributions to the physical observables come from a finite set 
of divergent Feynman integrals with the same number of internal momentum integrations 
(number of loops). Thus, one has to deal with very singular (ultraviolet and infrared) 
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Feynman integrals, and a correct mathematical apparatus suitable for calculational pur- 
poses is necessary. The dimensional regularization technique ('t Hooft and Veltman, 1972, 
1973; BoUini and Giambiagi, 1972; Ashmore, 1972; Cicuta and Montaldi, 1972) for ultravi- 
olet divergent Feynman integrals is based on the idea of integration over the space-time of 
noninteger dimension less than 4. In this case, the Feynman integrals become well defined, 
and divergences appear as poles in terms of the deviation from the physical space-time 
dimension of 4. The important property of dimensional regularization is that it pre- 
serves explicit gauge invariance and is very convenient for practical calculations. In fact, 
almost all recent progress in higher order analytical perturbative calculations has been 
made within the dimensional regularization, using 't Hooft's (1973) minimal subtraction 
prescription. 

The systematic study of strong interaction effects for the various observables in pro- 
cesses with large momentum transfer requires one to evaluate at least the first few coef- 
ficients in the perturbative expansion in terms of the strong coupling. Here the problem 
of calculating multiloop Feynman diagrams arises. The recursive type algorithm for an- 
alytical evaluation of one-, two- and three-loop massless, propagator type, dimensionally 
regularized Feynman diagrams has been given by Chetyrkin and Tkachov (1981), and 
Tkachov (1981, 1983a). This algorithm, together with the so called Gegenbauer x-space 
technique (Chetyrkin and Tkachov, 1979) allows one to evaluate an expansion in the Lau- 
rent series in £ = (4 — £))/2 of all massless propagator type Feynman diagrams up to 
the three-loop level, where D is the noninteger dimension of the space-time. The above 
algorithm is applicable to a wide class of problems up to the four-loop level. These are, 
for instance, calculation of renormalization constants, renormalization group functions, 
some of the cross sections and decay widths. We note that this algorithm deals only 
with propagator type massless diagrams. Nevertheless, due to the remarkable properties 
of dimensional regularization ('t Hooft and Veltman, 1972) and the minimal subtraction 
prescription ('t Hooft, 1973), namely, that the counterterms are polynomials in dimen- 
sional parameters within minimal subtraction (Collins, 1974; Speer, 1974; see also the 
textbook by Collins, 1984), a wide class of problems can be reduced to the evaluation 
of propagator type diagrams (Vladimirov, 1978, 1980). At high energies, in some cases, 
it is possible to neglect the masses of participating particles and consider the massless 
diagrams. The mass corrections of the type m^"/s", where s is the center-of-mass energy 
squared, can also be evaluated through the calculation of massless diagrams (see, e.g., 
Gorishny, Kataev and Larin, 1986; Surguladzc. 1989a. 1994a,b,c). Feynman graphs can 
also contain virtual heavy particle propagators regardless of the energy scale of the partic- 
ular process. If the masses of the virtual particles are much larger than the energy scale, 
one can neglect them, since their effcts are supressed by powers of large mass, according to 
the decoupling theorem (Appelquist and Carazzone, 1975). However, in some cases, such 
effects may not be entirely negligible (Soper and Surguladzc, 1994). The prescriptions 
to study asymptotic expansions of Feynman integrals in powers of w? / s can be obtained 
from Chetyrkin and Tkachov (1982), Tkachov (1983b,c, 1991, 1993), and Chetyrkin (1991; 
see also Smirnov 1990, 1991 and references therein). An exact general expression for one- 
loop, N-point, massive Feynman integrals has been obtained by Davydichev (1991), and 
Boos and Davydychev (1992). This expression contains the generalized hypergeometric 
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function and is complicated, except for some particular cases. An alternative method for 
massive Feynman integrals has been suggested by Kotikov (1991). 

In practice, the calculation of physical quantities within perturbation theory is very 
cumbersome and tedious already beyond the one-loop level, especially in realistic quntum 
field theory models, like QCD. However, the recursive type algorithms by Chetyrkin and 
Tkachov (1981) allow convenient implementation within algebraic programming systems 
like REDUCE (Hearn, 1973), SCHOONSCHIP (Veltman, 1967; Strubbe, 1978; Veltman, 
1989) and FORM (Vermaseren, 1989). Several computer programs were written in the 
last decade for analytical computation of multiloop Feynman diagrams. Among them we 
mention the programs which fully implement the above mentioned recursive algorithms. 
The program LOOPS (Surguladze and Tkachov, 1989a), written on the REDUCE system, 
calculates one- and two-loop massless, propagator type Feynman diagrams for arbitrary 
structure in the numerator of the integrand and for an arbitrary space-time dimension. 
The program MINCER (Gorishny, Larin, Surguladze and Tkachov, 1989), written on the 
SCHOONSCHIP system, and the program HEPLoops (Surguladze, 1992), written on the 
FORM system, calculates one-, two- and three-loop massless, propagator type diagrams. 
The status of the existing program packages has been discussed recently in Surguladze 
(1994d). The above methods, algorithms and computer programs allow one to make 
significant progress in high order analytical perturbative calculations of several important 
physical observables. 

The other outstanding problem in perturbative calculations is the renormalization 
group ambiguity of perturbation theory predictions. Indeed, starting from a certain order, 
the perturbative coefficients become scheme-scale dependent, while it is obvious that the 
calculated observable cannot depend on any subjective choice of nonphysical parameters. 
Several approaches have been suggested to deal with the scheme-scale ambiguity problem. 
Among them we consider the so called fastest apparent convergence approach (Grunberg, 
1980), suggesting one absorb the leading QCD corrections in the definition of the "effec- 
tive" running coupling. We will consider an approach based on the principle of minimal 
sensitivity of the physical observables to nonphysical parameters (Stevenson, 1981a,b), 
and Brodsky, Lepage and Mackenzie (1983) (BLM) method, suggesting one should fix the 
scale according to the size of the quark vacuum polarization effects. The commensurate 
scale relations by Brodsky and Lu (1994, 1995) allow one to make scale-fixed perturbative 
predictions without referring to the particular renormalization prescription. 

In the recent works, some authors try to predict the perturbative coefficients without 
calculating the relevant Feynman graphs. First, we mention the method by West (1991) 
which is based on the renormalizability, analyticity arguments, and the saddle point tech- 
nique. For comments on this work see Barclay and Maxwell (1992a), Brown and Yaffe 
(1992), Surguladze and Samuel (1992), and Duncan et al. (1993). The method of Samuel 
et al. (Samuel and Li, 1994a,b,c; Samuel, Li and Steinfelds, 1994a,b,c), based on Fade 
approximants, work surprisingly well for the large number of cases considered. However, 
a theoretical basis of this method is necessary. Recent developments have put Fade ap- 
proximant method on a more rigorous basis, which justifies its application to perturbation 
series in QED, QCD, Atomic physics, etc. This is discussed in recent papers (see, e.g.. 
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Ellis, Karliner and Samuel, 1995). An alternative method for estimation of higher or- 
der perturbative contributions can be obtained based on Stevenson's (1981a,b) approach 
(Surguladze and Samuel, 1993; Kataev and Starshenko, 1994). The important problem of 
large-order behavior of perturbation theory has been considered by Barclay and Maxwell 
(1992b), and Brown and Yaffe (1992). The same problem has been discussed during the 
past twenty years. The part of papers on the subject have been collected in the book 
edited by Le Guillou and Zinn- Justin (1990). The application of renormalon calculus in 
the study of the behaviour of perturbative QCD series is a subject of intensive discus- 
sions in the recent literature (see, e.g., Zakharov, 1992; Mueller, 1992; Lovett-Turner and 
Maxwell, 1994; Vainshtein and Zakharov, 1994, Soper and Surguladze, 1995). 

After a brief historical review, we turn to the main subject of the present article. 
Namely, we discuss the analytical high order perturbative calculations of several physical 
observables which have been completed recently with the help of the above mentioned 
methods, algorithms and computer programs. First, we consider the analytical calculation 
of R{s) in electron - positron annihilation at the four-loop level of perturbative QCD 
(Surguladze and Samuel, 1991a,b; Gorishny, Kataev and Larin, 1991), which turned out 
to be the most difficult among the problems of this type. This is the first and so far the 
only four-loop calculation of a physical quantity in QCD. |^ As a byproduct, the four-loop 
Rr in T decay (Gorishny, Kataev and Larin, 1991; Samuel and Surguladze, 1991) and four- 
loop QED (3 function (Surguladze, 1990; Gorishny, Kataev and Larin, 1990) have been 
evaluated. For earlier works, we mention, for instance, the calculation of the three-loop 
correction to R{s) in electron - positron annihilation (Chetyrkin, Kataev and Tkachov, 
1979; Dine and Sapirstein, 1979; Celmaster and Gonsalves, 1980), the calculation of the 
three-loop QCD f] function (Tarasov, Vladimirov and Zharkov, 1980) and the calculation 
of the three- loop anomalous dimensions of quark masses (Tarasov, 1982). We would also 
like to list some other three- and two-loop calculations. These are:the calculation of the 
total decay width of the neutral Higgs boson into hadrons at the three-loop level (Gorishny, 
Kataev, Larin and Surguladze, 1990, 1991b; Surguladze, 1994a,b), the calculation of the 
two- and three-loop Wilson coefficients in QCD sum rules (Surguladze and Tkachov, 1986, 
1988, 1989b, 1990), the calculation of the two-loop anomalous dimensions of the proton 
current (Pivovarov and Surguladze, 1991). So far only one five-loop calculation exists. This 
is the calculation of the five-loop renormalization group functions in 0^-theory (Kleinert 
et ai, 1991). 

The scope of the present paper is limited and we are not planning to review pertur- 
bative QCD. This has already been done and excellent reviews exist. We recommend, for 
instance, the recent work Handbook of Perturbative QCD by CTEQ collaboration (Brock, 
et ai), edited by G. Sterman (1993). Here we focus on a somewhat simplified descrip- 
tion of the key methods which allow one to perform analytical high order perturbative 
calculations up to and including the four-loop level. As an example, we will demonstrate 

^ This calculation was attempted earlier by Gorishny, Kataev and Larin (1988) but, unfortunately, 
errors were found. 

^ For the joint publication of the results of the two independent calculations of the four-loop QED 
/3-function, see Gorishny, Kataev, Larin and Surguladze (1991a). 



9 



the main points of the calculation of atot{e~^e~ — >■ hadrons), r(r~ — >• z/^ + hadrons), 
r{H hadrons) and the QED f3 function. We also outline the calculation of the Wilson 
coefficient functions of higher twist operators in the operator product expansion and dis- 
cuss various approaches to resolve the renormalization group ambiguity of perturbation 
theory predictions. 

The paper is organised as follows. In the 2nd section we introduce our notation and 
present some general relations. In this section we discuss the relevant methods and tools 
of perturbativc QCD. We briefly consider the necessary dispersion relation, the operator 
product expansion (OPE), the renormalization relations and the method for evaluation of 
the renormalization constants. We also discuss the main ideas of the method of projectors 
for calculating Wilson coefficients in OPE. In the 3rd section we evaluate the quantity 
r{H — > hadrons) at the three-loop level. In the 4th section we calculate the corrections 
to the correlation functions due to the nonvanishing quark masses. In the 5th section 
we describe the calculation of the Wilson coefficient functions of the dim=4 operators in 
the OPE of the two-point correlation function of quark currents. In the 6th section we 
describe the four-loop calculation of crtot(e'''e~ hadrons). Sections 7 and 8 are dedi- 
cated to the evaluation of r(r~ —>■ Ut + hadrons) and the QED j3 function respectively. 
In section 9 we discuss the problem of the renormalization group ambiguity of perturba- 
tive QCD results. As an example, we consider calculated quantities and use the known 
approaches to try to fix the scheme-scale parameter within the one parametric family of 
MS-type schemes. Next, we outline the original method of scheme-invariant analysis and 
optimization procedure by Stevenson (1981a,b). The paper ends with summarizing notes. 

II Calculational methods 

A Notation and general relations of perturbative QCD 

Throughout this paper we work within the standard model of strong interactions - QCD. 
For a review on QCD see, for example, Marciano and Pagels (1978), Mueller (1981), Reya 
(1981), and Altarelli (1982). For a textbook see, e.g., Yndurain (1983), Quigg (1986), 
Muta (1987), and Ellis and Stirling (1990). For the most recent source see, e.g.. Handbook 
of Perturbative QCD by CTEQ collaboration (Brock et al.), edited by G. Sterman (1993). 
The four-loop QED calculations will be discussed in section 8. 
The Lagrangian density of standard QCD is 

L{x) = -1/4(G«J2 - ^{d^Alf + Y^q^iid - mf)qf + g^^QfT^A'^Qf 

+a^c"^(a^5'^" + 5f/"*'=^yc^ (2.1) 

where G"^^ = d^A^ - d^^A^ + gf^^A^Al (a = 1, 2, 8) are the Yang-Mills field (Yang 
and Mills, 1954) strengths, A'^ and qj are gluon and quark fields, mj are the quark 
masses, are the Faddeev-Popov ghosts and is the gauge parameter. We use the 
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standard notation d = •y^d^ and A°- — 7'*^^. The index / enumerates the quark flavors, 
total number of which is N. The generators of the SUc(N) gauge group, the structure 
constants f"''"^ and d"'^'^ obey the following relations 



tab 



arpb 



tab 



(2.2) 



The eigenvalues of the Casimir operators for the adjoint [Na = 8) and the fundamental 
{Np — 3) representations of SUc(3) are 



3, Cj. = 4/3, and T = 1/2, rt'^'"^ = 40/3. 



(2.3) 



We use the standard QCD Feynman rules (see, e.g., Abers and Lee, 1973; Muta, 



1987). 



Propagators 

quark 



1 m+P f 



gluon 



1 5a6 
i P2 



ghost 



1 Sgh 

i P2 



Vertices 



quark-quark-gluon 



ghost-ghost-gluon 
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3-gluon 



{k-q)^, + g^,x{p-k)u] 



4-gluon 




The sum of all momenta coming in each vertex of the Feynman diagram is zero 
(momentum conservation). 

Factors 

(-1) for each closed fermion or ghost loop 

Statistical factors (for derivations see, e.g., 't Hooft and Veltman, 1973): 
I for each graph (subgraph) 



I for each graph (subgraph) etc. 
Integration 

Each loop corresponds to the integration / 

In general, the Feynman integral constructed according to the above rules is diver- 
gent. There are two kind of divergences. One, the so called ultraviolet (UV) divergence is 
due to large integration momenta and the other one - the so called infrared divergence is 
associated with the small integration momenta in the massless limit. The most convenient 
regularization of Feynman integrals is dimensional regularization ('t Hooft and Veltman, 
1972; BoUini and Giambiagi, 1972; Ashmore, 1972; Cicuta and Montaldi, 1972), where the 
space-time dimension is analytically continued from the physical value, 4, to a complex 
value D = A — 2e. In the limit £ — > 0, the divergences appear as poles 1/e, defining the 
counterterms. One of the remarkable properties of dimensional regularization is that the 
Ward identities implied by gauge invariance are maintained for arbitrary space-time di- 
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mension D, in contrast with the old Pauh-Villars regularization (Pauh and Villars, 1949). 
Another useful property is a convenience in practical multiloop calculations. Thus, in 
dimensional regularization we formally replace / J It is straightforward to 

extend the all necessary tensor algebra into D-dimensions. For example, g^^g^v = D, 
TrT/iTi/ = '^^^^gtiv, etc. For the complete list of formulae see, e.g., CoUins (1984) and also 
Narison (1982). Note, however, that the extension of the usual definition of the matrix 75 



is not straightforward. The totally antisymmetric tensor Scxp^iv is defined only in the 
four- dimensional space. In some cases the calculation of the quantities involving 75 is 
still possible within dimensional regularization. For a discussion of the problem of 75 
in dimensional regularization see Delbourgo and Akyeampong (1974), Trueman (1979), 
Bonneau (1980), Narison (1982), Collins (1984), and Larin (1993). For a calculation 
involving 75 within dimensional regularization see, e.g., Pivovarov and Surguladze (1991). 

In order to get finite physical quantities, the divergences in dimensionally regularized 
Feynman integrals, appearing as poles in l/s, need to be subtracted by adopting of some 
specific rule. This rule is usually called a renormalization scheme. Throughout this paper 
we use 't Hooft's minimal subtraction (MS) type scheme ('t Hooft, 1973). The subtrac- 
tion of divergences is equivalent to the redefinition (renormalization) of the parameters 
(coupling, mass and gauge fixing parameter) and fields in the original "bare" lagrangian 

= ii'^'Za.as, (i/V^TT = as) 

= mZ^, (2.4) 

// is a quantity of dimension of mass which is introduced within dimensional regularization 
in order to make an action dimensionless. Superscript "B" denote the unrenormalized 
quantity. We renormalize the gluon, quark and ghost fields analogously. Within the MS 
scheme the N-point Green function is renormalized in the following way 

r(pi, ...,pN,g,m,aG,ii) = Zrr^(pi, ...,pn, g^m^ac), (2.5) 

where is a polynomial in 1/e, and thus multiplying by Zr, we subtract only pole parts 
from the divergent F^. The evaluation of the renormalization constants Z will be discussed 
in the next subsections. 

It is easy to see that the fi parameter entered through the renormalization and hence 
the unrenormalized Green's function is independent of /i 

/x-^r^(pi, ...,pN,g,m,aa) = 0. 
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Using eq. ( p.5|) and expanding the full derivative we get the renormalization group equation 
in the following form 



d d d d 
l^^gj^ + Pias)as-g^ - 'jmi(ys)m— + PG{ois)-Q^ 7r 



r(Pi, ■■■,PN,'m,as,aG,fi) = 0. 



(2.6) 

The QCD renormalization group functions - the /5-function and the anomalous dimension 
functions - 7 are defined in the following way 



,daa 



dfx 



2 ' 



7m(as) = (2.7) 



7r(a.) 



/i^ dZr 



with bare coupling and mass fixed. In the present paper we use the renormalization 
group equation in the above form. The other forms are also known in the literature. The 
group properties of the renormalization was first discovered by Stueckelberg and Peterman 
(1953). The ultraviolet asymptotics of the Green function was studied by Gell-Mann and 
Low (1954) in quantum electrodynamics using the group of multiplicative renormaliza- 
tions. The renormalization group formalism was further developed in the original works 
by Bogolyubov and Shirkov (1955, 1956a,b). For the detailed monograph see Bogolyubov 
and Shirkov (1980). The renormalization group equation was studied by Ovsyannikov 
(1956), Callan (1970), and Symanzik (1970). For a recent historical review see Shirkov 
(1992) and references therein. 

The renormalization group /5-function and anomalous dimensions of quark masses are 
calculated up to the three-loop level (Tarasov, Vladimirov and Zharkov, 1980; Tarasov, 
1982). The QCD /5-function up to and including the three- loop level in MS type schemes 
is 

/?(«s) = -Po- - - H-f + 0{at), (2.8) 

TT TT TT 

where (Tarasov, Vladimirov and Zharkov, 1980) 



/3: 



2 
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The quark mass anomalous dimension up to and including three-loop level is 



Imias) = 70- + 7i(-)' + 72(-)' + 0{a% (2.9) 

TT TT TT 

where (Tarasov, 1982) 

70 = |Cf, 

71 = 1^ (ICf + f CfCa - f C^TTV) , 

72 = ^ [f^C^I - '-fClC^ + '-WCfCI - (46 - 48C(3))C|TiV 

-(^^ + 48C(3))Ci.CArA^ - ^CpT^N'' . 

As it was shown by Caswell and Wilczek (1974) and Banyai, Marculescu and Vescan 
(1974), the above renormalization group functions are gauge independent, which greatly 
simplifies their evaluation. In fact, the QCD /3-function and the quark mass anomalous 
dimension have been evaluated in the Feynman gauge = 1. We note that the per- 
turbative coefficients of the renormalization group functions arc the same within the one 
parametric family of the MS type schemes. Note also the independence of these pertur- 
bative coefficients on the quark masses by their definition within the MS type schemes. 



B Vacuum polarization function and Dispersion relation 

The vacuum polarization functions for various types of quark currents are crucial in the 
theoretical evaluation of total cross sections and decay widths. Indeed, for example, the 
quantity crtot(e"'"e~ hadrons), according to the well known optical theorem (see, e.g., 
the textbook by Bogolyubov and Shirkov, 1980), is proportional to the imaginary part of 
the function Il{—q^ + iO), defined from the hadronic vacuum polarization function 

^^^u{q) = i J e"^- < Tj^{x)UO) >o d^x = {g^^Q^ - g^g,)n(Q2)_l_. (2.10) 

Here, j^{x) = Qfqf'Jfj.qf, Qf is the electric charge of the quark of flavor / and = —q^ 
is the Euclidean momentum squared. The sum over all participating quark flavors is 
assumed in 11. The transverse form in the r.h.s. is conditioned by the conservation of 
electromagnetic currents. In this paper we also consider the two-point function of quark 
axial vector currents associated with the quantity V{Z hadrons) and two-point function 
of quark scalar currents associated with the quantity r(i7 — > hadrons) - the total decay 
width of the neutral Standard Model Higgs boson into hadrons. 

The renormalized vacuum polarization function obeys the dispersion relation 

4 /■°° R{ s] 

n(g2) = - / — -T^ds - subtractions, (2.11) 

6 J so S + 
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where 

D/ N (Ttot(e+e- ^ hadrons) 3 

Recall also that the muon pair production cross-section a{e^e^ — > fJ^^fJ^^) = Aira^/Ss, 
where a = e^/An is the electromagnetic fine structure constant. The above dispersion 
relation allows one to connect the experimentally measurable quantity R{s) to the Il{Q^) 
calculable perturbatively in the deep Euclidean region (Q^ is large compared to the typical 
hadron mass). For the discussion on theoretical calculability of R(s) see earlier references: 
Adler (1974), Appelquist and Politzcr (1975), De Rujula and Georgi (1976), Poggio, Quinn 
and Weinberg (1976), Shankar (1977), and Barnett, Dine and McLerran (1980). The 
combination of the idea of local duality in the dispersion relations (Logunov, Soloviov 
and Tavkhelidze, 1967) and the Operator Product Expansion technique (Wilson, 1969) 
became a basis of various versions of QCD sum rules (Shifman, Vainshtein and Zakharov, 
1979; Novikov et al., 1978, 1985; Krasnikov, Pivovarov and Tavkhehdze, 1983; Shifman, 
1992 and references therein). The methods of QCD sum rules are widely used to obtain 
quantitative information on the observed hadron spectrum and to extract the fundamental 
theoretical parameters. 

In practice, sometimes it is more convenient to introduce the Adler function (Adler, 
1974) 

Derivative here avoids an inconvenient extra subtraction in the r.h.s. 

The leading (parton) approximation of D{Q'^) in the zero quark mass limit coincides 
with R{s) 

D{Q') = 3J2Q}, (2.14) 
/ 

where the sum runs over all participating quark charges at the given energy. 3 stands 
for the number of different colors. The leading "non-QCD" contribution is completely 
free of ultraviolet divergences, while the Il{Q^) needs an additive renormalization even 
at the leading order. At higher orders of perturbative expansion of the Z)-function the 
ultraviolet divergences appear and one should employ a procedure (usually called renor- 
malization scheme) for their subtraction order-by-order. Because of ambiguity in the 
choice of subtraction scheme, the amplitude calculated within the perturbation theory 
depends on nonphysical parameters. Within the one-parametric family of the MS type 
schemes (t 'Hooft, 1973) such a parameter is usually called /i. Thus, up to power cor- 
rections, the £)-amplitude will be a function of log(/i^/(5^) and the strong coupling ag- 
On the other hand, since D is connected to the observable R{s), it can not depend on 
our subjective choice of nonphysical parameter fi. This can be achieved if the strong 
coupling becomes a function of providing independence of observables on the choice of 
parameter /i. Here, it is asumed that all orders of perturbation theory are summed up. 
Otherwise, if one considers a truncated series, the /i dependence remains. The problem of 
scheme-scale dependence and some possible solutions will be discussed later in this review. 
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The set of transformations which leave observables independent of renormahzation param- 
eters has a group character and forms the renormahzation group. The renormahzation 
group in renormahzable theories (hke QCD) fixes the dependence of the couphng on the 
yU-parameter. 

The function Di^Q"^) calculated in perturbative QCD within the MS type schemes 
obeys the renormahzation group equation 



D(/iVQ^m,a,) = 0. 



(2.15) 



Below we consider the limit of the massless light quarks and the infinitely large top mass 
which decouples (Appelquist and Carazzone, 1975). The solution of eq. ( |2.15| ) at /x^ = 
is 

D(/iVQ^ «,(/i)) = D(l, «,(Q)) = Has{Q)/7^y, (2.16) 



i>0 



where the as(/i^) is the running coupling, usually parametrized up to the three- loop level 
as follows 



as(/i^) 



1 



TT 



^ + 15^M log' L - Pi logL + P,Po - P!) + 0{L~'), (2.17) 



PoL PlL^ PIL^ 



where L = log(/i^/A^). Parametrization (|2.17| ) has the same form and the QCD /3-function 
coefficients are the same within the MS type schemes. The scale parameter A depends on 
the particular modification of the MS prescription. In fact, A is used to parametrize other 
versions of renormahzation prescription as well. It is shown by Celmaster and Gonsalves 
(1979) that the transformation relations valid to all orders between A's defined by any 
two renormahzation prescription can be deduced from a one-loop calculation. Comparing 
the bare coupling constants within different renormahzation prescriptions and using the 
results for the one-loop renormahzation constants and the property of asymptotic freedom, 
one obtains for, e.g., momentum subtraction (MOM) and MS schemes (Celmaster and 
Gonsalves, 1979) 



where 
A{ao,N) = CA 
+TN 
and the integral 



Ams — Amomcxp 



.-(,.-ln4.) + - + -/ 



2 4 
-(7e -ln47r) - - 



A{aa,N) 



(2.18) 



-ar 



I) - Y^«g(3 - /) 



1 s' 
— at 
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In X 







x + l 



-dx = 2.3439072.. 



(2.19) 



(2.20) 
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One note due to Stevenson (1981b, 1994) is in order. Despite its convenient form, the 
parametrization ( p.l7| ) produces an additional ambiguity due to the freedom with a partic- 
ular definition of A parameter, even when the renormalization prescription is already spec- 
ified. This problem was discussed by Abbot (1980), Shirkov (1980), Stevenson (1981b), 
Monsay and Rosenzweig (1981) and Radyushkin (1983). In fact, one can take advantage 
of this freedom in the choice of A and try to optimize the expansion in 1/L. Indeed, 
as was shown by Radyushkin (1983), if one takes 0.6A in eq. ( |2.17| ) instead of standard 
(Buras, Floratos, Ross and Sachrajda, 1977) A then the 1/L^ and terms contribute 
only a few percent for a reasonably wide range of /i. On the other hand, Stevenson (1981b, 

1994) has suggested to avoid the entire problem of ambiguity in the definition of A by 
abandoning the 1/L expansion and solving the renormalization group equation ( p.7|) for 
as and resulting transcendental equation numerically, using the truncated /3 function. 

According to the operator product expansion technique (Wilson, 1969), one can sep- 
arate perturbative and nonperturbative contributions to the function n(Q^). As shown 
by Shifman, Vainshtein and Zakharov (1979), this function can be represented in the 
following form 

n(<5^) = perturbation theory + ^ ^^^^ n2n " — ~ ~'~ i^istanton contributions, (2.21) 

n>2 ^ " 

where < 0„ >o denote vacuum condensates parametrizing the nonperturbative contri- 
butions and Cn{Q) are their coefficient functions. The last term in the above equation 
describes the instanton contributions, which, in the case of electromagnetic currents, was 
estimated to be small (Krasnikov and Tavkhelidze, 1982; Kartvelishvili and Margvelashvili, 

1995) . The coefficient functions of the condensates can be calculated within perturbation 
theory. High order perturbative corrections to the coefficient functions of dimension 4 and 
6 power terms have been calculated in Loladze, Surguladze and Tkachov (1984, 1985), 
Surguladze and Tkachov (1989b, 1990), Chetyrkin, Gorishny and Spiridonov (1985), and 
Lanin, Chetyrkin and Spiridonov (1986). In subsection E we discuss the method for evalu- 
ation of Wilson coefficient functions. Examples will be outlined in section 4. Note, that we 
consider the region of very high energies where, in fact, only perturbation theory contri- 
butions survive in n(Q^). The nonperturbative corrections could have some (small) effect 
in the case of, for instance, r lepton decay (see section 7). Note also that, in fact, the 
effects of neglected light quark masses are not entirely negligible in some phenomenological 
applications (see section 4). 



C Renormalization relations 

There are several approaches for the ultraviolet renormalization of Green's functions known 
in the literature. Throughout this paper we use 't Hooft's minimal subtruction method ('t 
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Hooft, 1971, 1973). For alternative prescriptions we refer to the works by Weinberg (1967), 
Gell-Mann and Low (1954), Callan (1970), Symanzik (1970), and Collins, Wilczek and Zee 
(1978). For an analysis of various renormalization methods see Collins and Macfarlane 
(1974). For a review see, e.g., Narison (1982), the textbook by Collins (1984) and references 
therein. We focus on the renormalization relations for the two point correlation function of 
quark currents relevant for the further evaluation of total cross sections and decay widths. 
It is known that the vacuum polarization function is renormalized additively 

n(/iVQ', «s) = n^(/iVQ', «.^) + Zn = finite. (2.22) 



The bare coupling is related to the renormalized one by the relation ( |2.4|) . The per- 
turbative expansion for Z^^ can be found based on eqs. ( p.7| ) and (pl8|), the MS definition 
of Za^ and the renormalization group equation 

/i^^«B ^ (2.23) 

We obtain 

. , _ , (?)^ (§-!)- (?)' (§ - ^ I) ^ (-^) 

In general, the polarization function depends on quark masses and we will need the relation 
between "bare" and renormalized masses up to 0{a1) (Tarasov, 1982) 

Within the minimal subtraction prescription ('t Hooft, 1973) the renormalization 
constant Z^ can be expressed as the following double sum 



Zn= E - (2-26) 

-i<fe<0 ^ ^ ^ 

Z>0 

where Zi^ are numbers. Furthermore, for the "bare" vacuum polarization function one 
has the following expansion in a perturbation series 

where the first index denotes the number of loops of the corresponding Feynman diagrams 
at the given order of ag. Substituting eqs. ( |2.27] ) and ( |2.22| ) into the definition ( |2.13| ), 
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after tlie renormalization of tlie coupling via ( |2.24 ) we obtain at = 



a. 



- ni,_i + — 



TT 



e 



+ { — 

TT 



- (3n3,-3 - 2/5on2,-2) 



+- (3n3,_2 - 2/?on2,-i) + (3n3,_i - 2/5on2,o) 



+ — 



3 r 



- 4n4,-4 - 6/?on3,-3 + 2p^U2,-2 



+- 4n4,_3-6/3on3,_2-/3in2, 



2,-1 



6/?on3,-i - f3iU2,-i + 2/32n2,o 



+- (4n4,_2 



+ (404,-1 - 6/3on3,o - /3in2,o + 2/5o'n2,] 



(2.28) 



Because of the renormahzation group invariance of D{fi'^/Q'^,as), in the above equation 
we take /x^ = to avoid unnecessary logarithms. The renormalized expression for the 
/^-function must be finite in the limit e 0. Thus the coefficients of pole terms must 
vanish identically. This implies relations between the perturbative coefficients of 11 and 
the QCD /5-function. First, we note that prior to any renormalization the leading poles 
must cancel at each order of in the sum of all relevant Feynman diagrams. As shown 
by the actual calculation, this happens in each gauge invariant set of diagrams. 

n4,-4 = n3,_3 = n2,_2 = o. (2.29) 

Moreover, from the cancellation of nonleading poles we get 

3n3,-2 - 2/3on2,_i = 0, 
4n4,„3 - 6/3on3,^2 + 2/5o'n2,_i = 0, (2.30) 

4n4,-2 - QPoUs,-i - (3iU2,-i + 2/?2n2,o = 0. 

The above relations provide powerful tests of the calculation at its intermediate stages 
and are crucial. 

From eq. ( |2.22| ) we see that fully renormalized Il{Q'^,as) must be finite. Thus, sub- 
stituting eqs. ( |2.24| ) - ( |2.27| ) and ( p.29| ) in eq. ( p.22| ) we obtain the following expression for 
the divergent part of n(/i^/(5^, «s) at /i^ = 



divnK) = -(ni,„i + Zi,_i) + ^ 

e IT 



-(n2,-i + Z2.-i] 

e 



4(n3,-2 - Po^2,-l + ^3,-2) + -(n3,-l - Po^2fi + Z3,-i] 

^(n4,-3 - 2/5on3,-2 + /3o'n2,-i + ^4,-3) 
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+4(n4,-2 - 2/5on3,-i + (3lii2,o - /3in2,-i/2 + Z4, 



+-(n4,_i - 2/?on3,o + (3lii2,i - /5in2,o/2 + Z4,„i) 



= 0. 



(2.31) 



The leading poles in are absent at each order of (^2,-2 = Z^ _^ = Z4 _4 = 0) except 
the zeroth order. Taking into account eq. (|2.3CI|) , we obtain the other set of relations 
between the perturbative coefficients of 11, Z and QCD /3-function 



32'3^_2 + — 0, 

2.^4,-3 + (^0^3 -2 = 0, 
4Z4,_2 + 2/3o^3,-l+/5l^2,-l 



(2.32) 



ni,-i — —^1,-1, 

n2,-l = —^2,-1; 
n3,_2 = —Z3^-2 — (3qZ2-1, 
^3,-1 = ~^3,-l + /^0n2,0, 

n4,-i = -Z4,_i + 2/5on3,o + /5in2,o/2 - /^^na,!, 

n4,-2 = -^4,-2 - 2/30^3,-1 - /3lZ2,-l/2 + P^n2fl, 
n4,-3 = —^4,-3 ~ 2/?o^3,-2 — /3o^2,-l- 



M,-3 



(2.33) 



In section 6, the above relations will be used in the calculations of the four- loop total 
cross-section in electron-positron annihilation. 



D Method for evaluation of renormalization constants 

We now discuss the evaluation of renormalization constants within 't Hooft's MS scheme 
('t Hooft, 1973), using Vladimirov's method (Vladimirov, 1978) and the so-called infrared 
rearrangement procedure (Vladimirov, 1980; Chetyrkin and Tkachov, 1982). 

To calculate the renormalizaton constant Z^ for the one-particle-irreducible Green's 
function F, it is convenient to use the following representation (Vladimirov, 1978) 

Zr = l- ICR'T. (2.34) 
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The operator /C picks out all singular terms from the Laurent series in e 



i i<0 

R' is defined by the recursive relation 

R'G = G — ^ ICR'Gi...lCR'Gn x G /(Giu...uGn}: 

(2.35) 

Gi 

where the sum runs over all sets of one-particle-irreducible divergent subgraphs Gi of the 
diagram G. G/(^Gi+...+Gn) is the diagram G with the subgraphs Gi, G„ shrunk to a point. 
In fact, R' is the ordinary Bogolyubov-Parasyuk /^-operation (Bogolyubov and Parasyuk, 
1955a,b, 1956, 1957; for a textbook see Bogolyubov and Shirkov, 1980) without the last 
(overall) subtraction. Thus, R' subtracts all "internal" divergences only and is connected 
to the ordinary i?-operation in the following way 

R={1- }C)R'. 

To calculate the renormalization constant Z in eq. ( p.22|) , one should write a diagram 
representation of 11 and apply ICR' to the corresponding graphs ( eq. ( |2.34| ) ) or, in other 



words, one should evaluate the counterterms for each graph. The benefit of using relation 
( p.34| ) is based on the fact that the ICR' for each diagram is a polynomial in dimensional 



parameters (Collins, 1974; Speer 1974). This fundamental property of the 't Hooft's 
minimal subtraction prescription is the basic idea of the various versions of the infrared 
rearrangement technique (Vladimirov, 1980; Chetyrkin and Tkachov, 1982). 

As an example, we demonstrate the application of the ICR' operation to the three- 
loop QCD diagram contributing to the 0{a1) total cross section for the process e+e~ 
hadrons. 
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R'l \^ -2ICR'l \ - 2ICR' 



KR'i \ = IC 



The benefit of using tlie ICR' operation besides its convenience in actual calculations is 
as follows. Using the fact that the result of K.R' operation is a polynomial in masses and ex- 
ternal momenta of the diagram, one can remove the dependence on the external momenta 
by differentiating (usually twice is sufficient) with respect to the external momentum and 
then setting the external momentum to zero. However, in this case infrared divergences 
appear. In order to prevent this, one can introduce a new fictitious external momentum 
as an infrared regulator flowing along some of the lines of the diagram (Chetyrkin and 
Tkachov, 1982). Alternatively, one can introduce a fictitious mass in one of the lines of 
the diagram as an infrared regulator (Vladimirov, 1980). An appropriate choice of the 
fictitious momentum can drastically simplify the topology of the given diagram. Both ver- 
sions of the so called infrared rearrangement procedure simplify the calculation and make 
it possible to evaluate counterterms to four- and five-loop diagrams. The main result of 
the application of the infrared rearrangement technique can be formulated as follows. The 
problem of calculating the counterterms of an arbitrary /-loop diagram with an arbitrary 
number of masses and external momenta within the MS prescription can be reduced to the 
problem of calculating some I — 1 -loop massless integrals to 0{s'^) with only one external 
momentum. In the later sections, the full calculational procedure will be demonstrated 
for a typical four- loop diagram contributing to the photon renormalization constant. 



E Evaluation of Wilson coefficient functions in operator product expansion 

In this subsection we briefiy discuss the problem of evaluation of higher twist operator con- 
tributions to the hadronic vacuum polarization function. Those contributions are relevant 
in the analysis of nonperturbative contributions in some processes (e.g., hadronic decay of 
the r-lepton). We use the Wilson operator product expansion technique (Wilson, 1969) - 
mathematical apparatus allowing a factorization of the short distance contributions, which 
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are calculable perturbatively and large distance effects which can be parametrized with 
the vacuum condensates ( Shifman, Vainshtein and Zakharov, 1979; Novilcov et ai, 1985). 
In the perturbative evaluation of Wilson coefficient functions, we rely on the so called 
method of projectors (Gorishny, Larin and Tkachov, 1983; Gorishny and Larin, 1987; see 
also Pivovarov and Tkachov 1988, 1993 and references therein). An actual calculation for 
the coefficient functions of the operators of dim = 4 has been done in the work by Lo- 
ladze, Surguladze and Tkachov (1984, 1985), and Surguladze and Tkachov (1989b, 1990). 
The present discussion is based mainly on those works. Below we demonstrate the above 
technique in the case of the coefficient functions of gluon and quark condensates. 

Consider the operator product expansion of the T-product of two quark currents in 
the deep Euclidean region, — = Q"^ ^ oo 

T{Q) = I J rf^xe*^^TJ(x) J(0) = Yl Q{Q)O,{0), (2.36) 

i 

where J are quark currents. Ci{Q) are c-number coefficient functions containing all depen- 
dence on Q. Oi are local operators forming in general a complete basis. If the currents J 
are gauge invariant then, after averaging over the vacuum, only gauge invariant operators 



contribute to the r.h.s. of eq. ( |2.36| ). However, the renormalization procedure mixes gauge 
invariant operators with non-invariant ones and one has to consider the complete basis of 
operators of the given dimension. The following set of operators of the dimension 4 

0, = {G;,Y, Oi = mfqfqf, Ol = Qfiid - mf + gT^A'^)qf, 

O4 = {d,r){d,c^) + (5^5"" + gr'^A';,)AlG% - gJ2qjT^A%, (2.37) 

/ 

ab I „ rabc ac\ b 



05 = 9^c"((9^r + ^7/""^^^ 



IC 



is closed under renormalization together with the "operator" ~ m"^ (Spiridonov, 1984; 
Loladze, Surguladze and Tkachov, 1984, 1985). Our aim is to calculate coefficient functions 
of gauge invariant operators Oi and O^- Note that ~ operators can be ignored because 
of the special structure of the renormalization matrix for the basis ( ^.37] ). The Feynman 
rules for the operators (|2.37|) are (Surguladze and Tkachov, 1990) 
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O4 2(5"^(pV"-p''p") (2.38) 

The operators of the basis ( |2.37| ) are renormahzed as follows 

= {ZohjOf, (2.39) 



where the superscript B marks the same operators as in (|2.37|) but built from the "bare" 



fields, masses and couplings. The structure of the renormalization matrix Zq has been 
studied by Spiridonov (1984). In the MS type schemes Zq has the following form (Surgu- 
ladze and Tkachov, 1990) 



(2.40) 
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where only the matrix elements A and B are relevant. 



A= 1 



e 



e \ e 
Inserting eq. ( |2.39| ) into the expansion ( |2.36| ) we get 



5=i^fi-^r\ (2.41) 



r{Q) = Y.C.{Q)Of{Zo)^,. (2.42) 



Following the method of projectors (Gorishny, Larin and Tkachov, 1983), we define the 
projectors vTj satisfying the orthogonality condition and vanishing on higher spin operators 

vraOf] = 5,„ (2.43) 

TTj [higher spin operators] = 0. 



Projectors tTj applied on the l.h.s. of eq. (|2.42| ) separate in the r.h.s. the coefficient functions 
we are interested in 

i 

We find the coefficient functions 

a(g) = E^.[^(^?)](^o')..- (2.45) 

3 

Our aim is to find the coefficient functions of gauge invariant operators Oi = {G^,y)'^ and 
O2 = mfq^qf. So, we need to construct the corresponding projectors tti and vTg. Let us 
represent TT^ clS 3i linear combinations of some "elementary" projectors Vj defined in the 
following way. 



p=m,f=0 



-pfp] = J—-^Tr( 

ANpdnif I Jp=mj.=o 

Hp] = i^i^Trrl } (2.46) 



ANp dp'' 



p=mf=0 



Na dp^ I J p=m;=0 
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where the parentheses contain the one-particle-irreducible Green function with one oper- 
ator insertion. In the case of Vi and P/ the traces are calculated over Lorentz spinor and 
color indices. 

Acting by the projectors Vj on the operators (|2.37|) we obtain 

Pi[Oi] = 8D{D - 1), Vi[0,] = AD{D - 1), p|[02^'] = 6ff,, 

P|[0{'] = -%, Vl[0(] = D5ff., V,[0,]=V,[0,]=2D, VrAO,] = D. (2.47) 

The results which are not shown in the above list are identicaly zero. From the definition 
( p.43| ) and eq. ( 2.47 ) we obtain the explicit form for the projectors vti and 



4 = H + ^Vi (2.48) 



Combining eqs. (p.40|) and ( p.41|) with eq. (|2.45|) we get our final expressions for the 



coefficient functions Ci{Q) and Cli^Q) (Surguladze and Tkachov, 1989,1990) 

Ci(g)=7ri[T(Q)](l-^), 

CiiQ) = 4[Tm - 7ri[T(g)]^^^. (2.49) 

The above expressions have a closed form and are valid at any order of perturbation theory. 
We note that T must be constructed with unrenormalized couplings and fields before one 
applies the projectors vTj. 

The general theory of Euclidean asymptotic expansions of Feynman integrals and 
the methods applicable to high order perturbative calculations have been developed in 
the works of Tkachov (1983b, 1983c, 1991, 1993), Chetyrkin and Tkachov (1982), and 
Chetyrkin (1991) (see also Smirnov 1990, 1991 and references therein). The technique 
developed in these works allows one to derive operator product expansions in the MS- 
scheme for any Feynman integral. For more general discussion and further details we refer 
to the above works and also to the original calculations (Surguladze and Tkachov, 1989, 
1990). In section 5 we present a short description of the calculation of the coefficient 
functions of gluon and quark condensates up to 0{al). 
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Ill r{H hadrons) to O^a]) 



A The decay rate in terms of running parameters 

In this subsection, using the above methods we calculate the 0(0;^) corrections to the 
total hadronic decay width of the Standard Model Higgs boson in the massless quark 
hmit (Gorishny, Kataev, Larin and Surguladze, 1990, 1991b; Surguladze, 1994b). 



FIG. 1. The process H — > hadrons 

The standard SU(2)xU(l) Lagrangian density of fermion- Higgs interaction is 

L = -gyqfQfH = -{V2GFy/^mfqfqfH = -(v^G^ (3.1) 

The decay width of a scalar Higgs boson to the quark- ant iquark pair is determined by the 
imaginary part of the two-point correlation function 

n(g2 ^-s,mf)^ij e*«^ < Tjf{x)jf{0) >o d^x (3.2) 

of the quark scalar currents jf = uifq^qf in the following way 



^H^qfqj = ^^^ ImH(g + iO, m/) 



(3.3) 



Mh is the Higgs mass. The total decay width will be the sum over all participating 
(depending on Mh) quark flavors 

r(i7 hadrons) = ^ Tn^qf-qf (3.4) 
f=u,d,s,... 
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We follow the work by Gorishny, Kataev, Larin and Surguladze (1990) and in analogy 
to the vector channel introduce the Adler function (Adler, 1974) 

The derivative avoids the additive renormalization of 11. In fact, it is possible to proceed 
without the introduction of D-function and deal directly with the correlation function 11 
(Surguladze, 1994b). Indeed, we are interested in Imn(s + iO,mf). Since the overall MS 
renormalization constant has no terms like [log n"^ / Q'^)"' / e'^ , its imaginary part vanishes 
identically. The abscence of the pole logarithms in renormalization constants is a general 
feature of MS type schemes. 

The Z)-function obeys the homogeneous renormalization group equation 

( o d d d \ 

+ /3(a^)a^- lm{ois)^ \ D{i/ IQ'^.mf.as) = (3.6) 
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The QCD /3-function and the mass anomalous dimension 7^ are known up to the three 
loop approximation and have been given in the previous section. The plan for evaluation 
of Tn^gf-qj, is as follows. First, we write the diagram representation for n((5^,m/) ac- 
cording to the standard Feynman rules up to the desired loop-level. Second, we evaluate 
the Feynman diagrams using the dimensional regularization and renormalize the coupling 
and quark masses within the MS renormalization prescription. Finally, to get the decay 
rate, we analytically continue the result for the D-function obtained from eq. (|3.5|) from 
Euclidean to Minkowski space. Following the above plan, we now demonstrate the calcu- 
lation of Tn^qj-qj up to the 3-loop level. First of all, note that the correlation function 
n and the related D-function depend on quark masses. The algorithms for evaluation of 
the 3-loop Feynman diagrams constructed with the propagators of massive particles has 
not yet been developed. However, in the deep Euclidean region (Q^ — > 00) it is possible 
to simplify the calculation using the expansion in terms of the small parameter m^/Q^ 

.2 . 



1 / TTi 

^(g^m;) = ^(g2) + o -^). (3.7) 
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Such an expansion is legitimate since we consider a Higgs boson much heavier than the 
typical hadronic mass scale. In this section we calculate the first term in the above 
expansion and the related decay rate. This is equivalent to the assumption that all five 
quarks are massless and the top quark decouples {rrit — > oo). 

The diagrammatic representation for 11 in somewhat symbolic form looks like 
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Next, we evaluate one-, two- and three-loop massless Feynman diagrams. By simple power 
counting, it is easy to find that in general the above diagrams are UV divergent. The 
unrenormalized contribution from a typical three-loop diagram in the MS renormalization 
scheme (Bardeen, Buras, Duke and Muta, 1978) reads 
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where is the /-flavor quark mass originating from the quark mass dependence of the 
Yukawa coupling. C(3), C(4) and C(5) are ordinary Riemann ^-functions. The number 2 
in front of the diagram stands for the symmetry factor. The algorithms for the evaluation 
of propagator type one-, two- and three-loop massless Feynman diagrams have been given 
by Tkachov (1981, 1983a) and Chetyrkin and Tkachov (1981). For the description of the 
algorithms see also Gorishny, Larin, Surguladze and Tkachov, 1989. The results given 
in this section were reobtained with the help of the program HEPLoops (Surguladze, 
1992) and the previous results (Gorishny, Kataev, Larin and Surguladze, 1990, 1991) were 
independently confirmed (Surguladze, 1994b). 

As one can see, each three-loop diagram in general may contain a pole with power 
< 3. In the vector channel, after summing the results for all diagrams with an appropriate 
symmetry and SU(N) group factor, the leading pole cancels. This is the consequence 
of the conservation of electromagnetic currents. In the scalar channel, the leading poles 
remain in 11. This is related to the quark mass dependence of the coupling. 

Evaluating the unrenormalized correlation function (|3.2| ) and using the definition 
(B.5|), we obtain the unrenormalized D-function in the massless limit. 
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+ 0(a^)}. (3.9) 



The above expression requires the renormahzation of the strong couphng ( eq. ( 2.24 



and the muhiphcative renormahzation ( eq. ( p.25| ) ) originating from the quark mass 
dependence of the Yukawa couphng. 

Expanding the factors {^j^/Q'^Y^ in terms of e and performing the renormahzations 
of the couphng and the quark mass, we get a finite analytical expression for the D-function 
in the MS scheme 
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For standard QCD with the color SUc(3) symmetry group, the analytical result for the 
D-function reads (Surguladze, 1989d) 
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This completes the evaluation of the correlation function of the two scalar quark currents 
in the massless limit at the three-loop approximation. 

There is one crucial test of this calculation based on the renormahzation group con- 
straints. The solution of the renormahzation group equation (|3.6| ) can be conveniently 
rewritten as follows 
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Applying the differential operator fj,'^d/dfi'^ to both sides of eq. (|3.12|) , taking into account 
the renormahzation group invariance of the D-function and eqs. ( p.8|) and (|2.9| ), we obtain 
to 0{as) 

an = 27oaoo, (3.13) 
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to 0(a2) 

aai = 271 aoo + (/^o + 270) Oio, 
a22 = {l3o + 270)^ = (/5o + 27o)7o«oo, (3.14) 

and to 0{al) 

031 = 2(/?o + 70)020 + (A + 27i)aio + 272000, 
032 = (/5o + 70)021 + (A + 271)^ = (/5o + 7o) [271^00 + (/5o + 27o)aio] + (/3i + 271)70000, 

2 2 

033 = g(/5o + 7o)a22 = -7o(/5o + 7o)(/3o + 270)000- (3.15) 

The relations ( ^.13| ) and ( |3.14| ) provide a powerful check of our calculation, while the re- 
lations ( ^.15| ) allow one to evaluate the log terms to 0(q;^), without explicit calculations 
of the corresponding four-loop diagrams. With those relations, the information available 
at present, namely the QCD /3-function, mass anomalous dimension and the two-point 
correlation function up to the three- loop level is fully exploited. In fact, similar relations 
can be derived for the correlation function 11. However, the renormalization group equa- 
tion for n is not a homogeneous one and the anomalous dimension function up to the 
corresponding order of is necessary. 

We evaluate the decay rate of the neutral Higgs boson into a quark antiquark pair by 
analytical continuation of D (/i^/Q^, ""^/(/x), as(/^)) from Euclidean to Minkowski space. 
The total decay rate can be obtained by summing up over all participating quark flavors. 
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The Riemann function ^(2) = n'^/Q arose from the analytical continuation of the log^ /^Its/ 
term and ({3) = 1.202056903. The procedure of analytical continuation and the appear- 
ance of invariant additional contributions have been discussed in several earlier works 
(Krasnikov and Pivovarov, 1982; Pennington and Ross, 1982; Radyushkin, 1982; Pivo- 
varov, 1992a). Note that in some cases those additional corrections are large and affect 
the result significantly. This is especially true for the total cross section in the process 
e+e" — >■ hadrons. To minimize such corrections it was proposed, for instance, to redefine 
the expansion parameter (Pennington and Ross, 1982; Radyushkin, 1982). 



B The decay rate in terms of pole quark mass 

For the heavy flavor decay mode of the Higgs, it is relevant to parametrize the decay 
rate in terms of quark pole mass (see, e.g., Kniehl, 1994a). Let us rewrite the result for 
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VH_,qj-q^ in terms of pole quark mass, assuming that heavy quark is not exactly on-shell. 
This subsection is based mainly on recent findings (Surguladze, 1994a,b). 

Solving the renormalization group equation for the quark mass - eq. ( p.7|) , we obtain 
the following scaling law for the running quark mass: 
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In the above equation all appropriate quantities are evaluated for active quark fiavors. 

can be determined according to the scale of M^- At present one usually considers 
N = 5. 

For the running coupling we obtain the following evolution equation to O(a^) (Sur- 
guladze, 1994b) 
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where the superscript (n) indicates that the corresponding quantity is evaluated for 
(n) numbers of participating quark flavors. Conventionally (see, e.g., Marciano, 1984) A^ 
(n) is specified to be the number of quark fiavors with mass < /i2 (< /Ui). However, eq. 
( ^.19| ) is relevant for any n < N and arbitrary /xi and fi2, regardless of the conventional 
specification of the number of quark fiavors. The log mi / fii terms are due to the "quark 
treshold" crossing effects and the constant coefficients 1/6 = Pq'' —Pq^\ 19/24 = — 

(k) 

Pi represent the contributions of the quark loop in the /3-function. The sum runs over 
N — n quark fiavors (e.g., / = 6 if n = 4 and A^ = 5). Note that mi is the pole mass of the 
quark with fiavor /. For the on-shell definition of the quark masses eq. ( |3.19D changes - the 
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constant —25/72 should be substituted by +7/12. The above equation is derived based 
on eq. ( p.l7| ), the QCD matching conditions for at "quark thresholds" (Bernreuter and 
Wetzel, 1982; Marciano, 1984; Barnett, Haber and Soper, 1988; Rodrigo and Santamaria, 
1993) and the one-loop relation between on-shell and pole quark masses. Eq. ( |3.19|) is 
consistent with the QCD matching relation at m/(m/) (Bernreuter and Wetzel, 1982) 

aff-'\mf{mf)) = /)(m;(m;)) + (af /)(mKm;)))3(CV9 - 17Cj./96)y (3.20) 

Here and below Nf is the number of quark flavors u,d, f. Note that the nonlogarithmic 
constant at 0{al) in eq. ( |3.19|) will not contribute in further analysis. 

Next, using the scaling properties of the MS running mass and eq. ( p. 191) , one obtains 
the following matching condition 
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where the constant terms are: 1/12 = 7o(/?o'' — /9o'^'')/2, 5/36 = 7i ' — 7i 
2/9 = CpipH"^^ - pH"^). In general, the 6{mf, mf) is the finite contribution of the single 
virtual heavier quark with mass rrifi, entering when one increases the number of flavors 
from — 1 to (one can also consider the particular case mf> = ruf). 

From the two-loop on-shell quark mass renormalization one has (Broadhurst, Gray 
and Schilcher, 1991) 
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L±{r) can be evaluated for different quark mass ratios r numerically. We relate the MS 
quark mass mj(mj) to the pole mass rrif using the O(a^) on-shell results of Broadhurst, 
Gray and Schilcher (1991) 
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The first four terms in Kj represent tlie QCD contribution witli Nj massless quarks, wliile 
tlie sum is tlie correction due to tlie Nj nonvanisliing quark masses. 

Combining eqs. (|3.17| ), (|3.18|) and eqs. ( p.l9| )-( ^l2^ ), one obtains tlie relation between 
tlie MS quark mass nif^MH) renormalized at Mh and evaluated for the A^-flavor theory 
and the pole quark mass rrif (Surguladze, 1994b) 
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Note that is specified according to the size of Mh and has no correlation with the 
quark mass m/. Thus, for instance, one can apply eq. ( p.26| ) to the charm mass m^\MH) 
evaluated for five-flavor theory. 

Substituting eqs. (|3.26|) , ( p.25| ) and appropriate /^-function and mass anomalous 
dimension coefficients (see section 2) into eq. (|3.16 ), one obtains the decay rate in terms 
of the pole quark masses 
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Recall, that at the beginning we have neglected terms which are suppressed by pow- 
ers m'j/Ml. Such corrections to the decay rate, in general, may not be entirely negligible 
and have to be taken into account in precise numerical analyses. Presently those correc- 
tions due to the nonvanishing quark masses have also been calculated. For the explicit 
results, we refer to the original works (Surguladze, 1994a,b; Kniehl, 1995a; Chetyrkin and 
Kwiatkowski, 1995). In the next section we give the results for the quark mass corrections 
to the correlation functions 11. 

The full analytical result for the decay rate of if g/g^ in terms of pole quark masses, 
including the leading order (two-loop) QCD corrections has been obtained independently 
by several groups: Braaten and Leveille (1980), Inami and Kubota (1981), and Dreess and 
Hikasa (1990). In the work by Sakai (1980) the two-loop result has been obtained in the 
zero quark mass limit. 
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and the Spence function is defined as usual 
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The expansion of the r.h.s of eq. (|3.28|) in a power series in terms of small m'j/Ml has the 
following form 
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where the periods cover higher order terms ~ {mj / MhY^ ■, k = 2, 3... One can see that the 
leading terms agree with the result ( |3.27| ). 

Numerically the MS high order QCD corrections for the considered process are large 
and reduce the decay rates by about 40%. 
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IV Quark mass corrections to the correlation functions 

In the previous section we neglected all quark masses in the corresponding Feynman 
diagrams in comparison with the momentum scale of the problem. In other words, we 
have calculated the leading term in the expansion in terms of small m^/s (for the Higgs 
boson decay, s = Ml) in the limit of infinitely heavy top quark, 00. However, in 

the real world quarks are massive and the leading term in the above expansion may not 
always give a satisfactory approximation. On the other hand, starting at O(a^), virtual 
heavy quark can also appear in certain topological type of Feynman diagrams (Fig. 2 and 
Fig. 3) regardless of the momentum scale of the problem. 
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FIG. 2. O(a^) Feynman diagrams responsible for the virtual heavy quark contribution 



FIG. 3. O(a^) Feynman diagrams responsible for the contribution due to the 

top-bottom mass splitting 

According to the decoupling theorem (Appelquist and Carazzone, 1975), virtual 
quarks much heavier than the momentum scale of the problem decouple. However, for 
instance, in the process of Z boson decay the effect of the top quark may not be entirely 
negligible since rrit is not much greater than Mz- A similar role could be played by the 
charm quark in the hadronic decay of the tau-lepton. The evaluation of the virtual top 
quark contribution (Fig. 2) to the decay rate Z —>■ hadrons and related quantities has 
been done in Kniehl (1990), Soper and Surguladze (1994), and Hoang, Jezabek, Kiihn and 
Teubner (1994) without using large or small mass approximations. The correction turned 
out to be moderate and in good agreement with the results obtained with the help of the 
large mass expansion technique (Chetyrkin, 1993a). The contribution of the diagrams in 
Fig. 2, in the presence of a virtual heavy quark, to the two-point correlation function of the 
electromagnetic quark currents has been evaluated previously by Wetzel and Bernreuther 
(1981). Kniehl and Kiihn (1989, 1990) have calculated the O(a^) correction to the decay 
rate Z hadrons due to the large mass splitting in the top-bottom doublet (Fig. 3). This 
correction turned out to be large and important. 

In this section we consider only the leading correction in the expansion in terms 
of small quark mass. For the calculations of virtual heavy quark contributions we refer 
the reader to the above mentioned original works (see also Kniehl, 1994b, 1995b). The 
discussion in this section is based on the works by Surguladze (1994a,b,c). 

Let us expand the full two-point correlation function, defined by eq. ( |2.10| ) in the 



vector channel and by eq. ( ^.21 ) in the scalar and pseudoscalar channels, in powers of 



m'l / in the "deep" Euclidean region 



^ J U{Q',mf,m,) =U,{Q') + -Lu^.iQ') + Yl ^n^^(Q2) + (4.I) 

^ v=u,d,s,c,b ^ 

where d = in the vector channel and d = 1 in the scalar and pseudoscalar channels. 
The last term in the above expansion is due to the Feynman diagrams containing a vir- 
tual fermionic loop. Note however that in the vector channel the contribution from the 
diagrams in Fig. 3 vanishes according to Furry's theorem (Furry, 1937). 
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In order to evaluate the coefficient functions in the r.h.s of eq. (|4.1| ), it is sufficient 
to write the diagrammatic representation for Il{Q^ , , m^) up to the desired level of 
perturbation theory and apply the appropriate projector. To O(a^) one has 
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where n, /c = 0,l,n + /c<l, and superscript "B" denotes the bare quantities. The mass 
renormalization constant Zm = m^/mf can be obtained from eq. (|2.25|) . The Feynman 
diagrams contributing to the 11^2 the same as for the calculation of Hi (see eq. 

( p.8|) ) but with massive fermion propagators. The calculations of all one-, two- and three- 
loop diagrams have been done using the program HEPLoops (Surguladze, 1992). 

The obtained expressions for Ilj at each order of are polynomials with respect 
to 1/e and log/i^/Q^. The poles can be removed by an additive renormalization. We 
note that there are no terms like {l/e"')(\og^^/Q'^)'^. They appear only at higher orders 
~ m'jm'j/Q'^ and represent infrared mass logarithms. The corresponding prescription 
similar to the Bogolyubov ultraviolet i?-operation has been worked out in the work by 
Chetyrkin, Gorishny and Tkachov (1982), Tkachov (1983b, c), and Gorishny, Larin and 
Tkachov (1983). (see also Tkachov, 1991, 1993 and references therein). The infrared 
mass singularities have been studied earlier by Marciano (1975). In the present paper we 
consider only the terms ~ m'j / which are sufficient for most of the phenomenologically 
interesting applications. 

In the vector channel we obtain the following MS analytical result (Gorishny, Kataev 
and Larin, 1986; Surguladze, 1994c) 
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The contribution to the physical process, in particular to the decay rate of Z 
hadrons can be obtained simply by taking the imaginary part in the r.h.s. of eqs. 
and (O) at Q"^ = -s + iO. We note, that the U 



2 and turned out to be finite. No 
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overall subtraction is necessary. Moreover, one can see that the imaginary part or the 
contribution to the decay rate vanishes at the parton level. This can be checked by the 
calculation of the parton contribution in the vector channel with explicit dependence on 
quark mass. Indeed, calculating the trivial fermionic loop we obtain 
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Taking the discontinuity under the integral and then evaluating the trivial integral with 
the G function, we obtain 
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The 
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1 1 contribution to the Adler D-function can be obtained from eqs. 



and 



^ by differentiating with respect to . 

There is some confusion concerning the above results in the literature. Initially, the 
corrections ~ m)/Q'^ in the vector channel have been calculated by Gorishny, Kataev 
and Larin (1986). Later, in the similar calculations (Surguladze, 1989a), a slightly differ- 
ent result was obtained, which was confirmed in further publications (see, e.g., Kataev, 
1990, 1991). However, in the recent works (Chetyrkin and Kwiatkowski, 1993; Surguladze, 
1994c), the initial result of Gorishny, Kataev and Larin (1986) has been confirmed. Un- 
fortunately, in the analysis of the mass corrections to the Z decay rates (Chetyrkin and 
Kiihn, 1990) the incorrect result was used. Fortunately, the main conclusions of Chetyrkin 
and Kiihn (1990) are not affected. Summarizing, we note that the results (|4.3| ) and (^!^ ) 
(Gorishny, Kataev and Larin, 1986; Chetyrkin and Kwiatkowski, 1993; Surguladze, 1994c) 
seem now to be reliable. 

In the scalar channel the result for the standard SUc(3) gauge group reads (Surgu- 
ladze, 1994b) 



n 2 ( ^ a 



+ 



+ - 

V vr 



TT 



1 



4vr2 



12 + 9 log 



2 2 

94 - 36C(3) + 60 log ^ + 18 log' ^ 



17245 
16 



15E2c(3)-3C(4) + 5?5c(5) 



.(T-.«c(3,),o4.iH!vf.| 



-N 



817 
"2? 



6C(3) + 



/313 
\l2 " 



6C(3))log^ + ^log2^ + log^ 



^3 f^jAS 



Q 

+ "simple poles" 



(4.7) 



39 



47r2 



TT 



8 u— 

- + 6 log ^ + "simple pole" 
3 



(4.8) 



where under the "simple pole" we mean number/e*^ with no dependence on log/i^/Q^. 
The "simple poles" have no imaginary part and consequently will not contribute to the 
observable quantities at the given order of a^. Note that the 11^2^ in eq. (^]8|) does not 
include the contribution from the triangle anomaly type graphs pictured in Fig. 3. Those 
graphs make the following additional contribution to 11 in eq. ( [4.1|) (Surguladze, 1994b) 
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The above results are relevant for the decay rate of the standard model Higgs boson 
into a quark antiquark pair, calculated in the previous section in the massless quark 
limit. Corrections ~ m'j/Mfj can be obtained from eqs. ( [4.7|) , ( [4.8|) and (^4.9|) (Surguladze, 
1994b). 

In the pseudoscalar channel we define the quark currents as jf = rrifqji'j^qf. We 
also define the 75 matrix in D-dimensional space-time as an object with the following 
properties 

{75,7/.} = 0, 7575 = 1. (4.10) 

The above definition causes no problems in dimensional regularization when there are two 
75 matrices in a closed fermionic loop. We obtain (Surguladze, 1994a) 
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The result for the pseudoscalar channel is relevant, for instance, for the decay rates of 
the minimal supersymmetric version of the Higgs particle into a quark antiquark pair (see 
Surguladze, 1994a). 



Finally, we present the results of calculation of the 



m 



c/Q' corrections to the 



correlation function in the axial channel (Soper and Surguladze, 1994; Surguladze, 1994c). 
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We use the following definition of the correlation function 



I J d'xe"^^ < TjI{x)£{0) >o= g,.Q'U{Q,mf) - Q^Q,U'{Q,mf), (4.13) 

where = Qf'jfj.'y^qf- Note that in the axial channel the correlation function is not 
transverse in contrast to the vector channel. However, for the decay rate of the Z-boson 
only the ~ (7^1/ part in eq. (|4.13| ) is relevant. 



The expansions of 11 and 11' in terms of small m'j/Q'^ has the same form as in the 
vector channel ( eq. ( [4.1| ) ). The coefficient functions in this expansion can be calculated 
according to eq. (|4.2| ) in the vector channel. In the calculations of one-, two- and three- 
loop Feynman diagrams the program HEPLoops (Surguladze, 1992) was used. The final 
results for the SUc(3) gauge group read (Soper and Surguladze, 1994; Surguladze, 1994c) 
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The results given in this section can be tested using the renormalization group. 
Namely, the relations similar to eqs. ( ^.13| ), ( ^.14] ) and ( pj.lS] ) can be obtained here (Sur- 
guladze, 1994a,b,c). In fact, in the vector channel, one can obtain the O(a^) logarithmic 
terms without actual calculation of the corresponding four-loop diagrams. On the other 
hand, the leading logarithmic terms in H-function form the corresponding contribution to 
the decay rates of, for instance, the Z-boson (Chetyrkin and Kiihn, 1990; Chetyrkin, Kiihn 
and Kwiatkowski, 1992; Surguladze, 1994c). In the axial channel the situation is more 
complicated. Here, because the renormalization group equation similar to eq. ( |3.6|) is no 
longer a homogeneous one, the renormalization group approach is restricted to 0{al). 
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V Two-loop coefficient functions of dim = 4 power corrections 



In this section we outhne the calculations of the two-loop coefficient functions of dim = 4 
power corrections. We consider the contributions which appear in the short distance 
expansion of the correlation function of two flavor- diagonal vector, scalar and pseudoscalar 
currents constructed from light quark fields. The methods and corresponding references 
are given in the earlier sections. The corrections for the vector channel have been evaluated 
in Loladze, Surguladze and Tkachov (1984, 1985) and Surguladze and Tkachov (1988). In 
the scalar and pseudoscalar channels, the calculation has been done in Surguladze and 
Tkachov (1990). The calculation for vector and axial vector channels has been done in 
Chetyrkin, Gorishny and Spiridonov (1985), where the previous results for the vector 
channel have been confirmed and the calculation was extended for flavor non-diagonal 
currents as well. The three-loop correction to the coefficient function of gluon condensate 
in the scalar channel has also been computed in Surguladze and Tkachov (1989b). For the 
calculation of dimension 8 terms in the operator product expansion see also Broadhurst 
and Generalis (1985). Here we follow the work by Surguladze and Tkachov (1990). 
Consider ffist the T-product of flavor diagonal vector currents of light quarks 



riliQ) = I J d'xe''^''Tjj:'{x)Jl'iO), 



(5.1) 



where J^' = Qf^y^qf- Taking into account the current conservation and operator product 
expansion technique (Wilson, 1969) for large momentum transfer (Q^ — * oo) we write 



^/JW) = (^?..Q'-Q.Q.){c^o + ^[c^G<g')(G^J' + E^^£(g')m/,g;,g/, 



+ ■ 



(5.2) 



where Cq is the coefficient function of the unity operator including the terms ~ mj/Q"^ 
discussed in the previous section. The period covers the operators of higher twists. For 
the scalar and pseudoscalar channels the transverse factor in the above equation is absent. 
To simplify the calculation, we contract over the Lorentz indices fi and u. Then the 
expressions for Cj defined in eq. ( |5^ ) coincide with the ones in eq. ( p.49| ) if T{Q) is 
replaced by T^^{Q'^)/{D — 1)(5^, where D = A — 2£. Let us rewrite eqs. ( |2.49| ) in a 
somewhat symbolic diagrammatic representation to 0{a1). 
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The total number of two-loop graphs contributing to Cc2 is 30 and to C^^ is 38. 
There is a simple rule for generating the appropriate graphs at 0(a"). One should take 
the graphs contributing to O the unity operator and disconnect one fermion 

line in all possible ways for the coefficient function C^^. For the coefficient function C^a 
it is necessary to write all the diagrams with one disconnected gluon line (relevant for 
the projector Pi), all the diagrams with one disconnected ghost line (relevant for the 
projector P4) and all the diagrams with disconnected gluon-ghost-ghost vertex (relevant 
for the projector V^). To see this, recall eqs. ( |2.46[ ). Acting with the projectors Vj on 



the appropriate diagrams, the calculations are reduced to the evaluation of one- and two- 
loop propagator type massless Feynman integrals. In the original calculation (Loladze, 
Surguladze and Tkachov, 1984, 1985; Surguladze and Tkachov, 1988, 1990) all Feynman 
integrals have been evaluated analytically using the REDUCE (Hearn, 1973) program 
LOOPS (Surguladze and Tkachov, 1989a). 

The MS results for the projectors Vj in the vector channel read 



+ 0{al) \ (5.4) 



nril] = ^C^^^{48-325 + ^[cH-12) + C^(^-42 + 72C(3)) 

^4[T/;] = ^,Cp^Ca[^)\1 - 9 + 12C(3)) + 0{al), (5.5) 

n[r/;]=0 + O(a^), (5.6) 

(Vi^f + Ip/^/') [r,/'] = ^,CpT[^)\^ - 60 + 96C(3)) + Oia% (5.7) 
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The vanishing of V^lTj^] at the two-loop level is the consequence of gauge invariance, 
as was shown by Spiridonov (1987). Combining eqs. ( p.48| ) and ( p.49| ) with the above 



results and renormalizing the bare coupling via eq. ( |2.24| ) we obtain MS 0(0;^) analytical 



expressions for the coefficient functions in the vector channel (Surguladze and Tkachov, 
1990) 
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Cifm = ^(v)'^^^(-^ + 2C(3)) +0(aD. (5.11) 

The above results are gauge invariant. This statement was checked by straightforward 
calculation in an arbitrary covariant gauge up to the term ~ e (Surguladze and Tkachov, 
1990). The dependence on the gauge parameter canceled. Thus, it is simplest to do the cal- 
culation in the Feynman gauge. For simplicity, we have omitted the terms ~ \og{fi^/ Q"^) , 
taking = Q^. The dependence on yU can be restored via the renormalization group 
(see below). Note that the coefficient function C^q is due to the diagrams pictured in 
Fig. 4 with disconnected fermion lines of the virtual loop (see also Fig. 2). 



FIG. 4. Two-loop diagrams forming C^^^ 



Specifically for QCD with the SUc(3) symmetry group we obtain 
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(5.13) 
(5.14) 



Note the very large 0{a1) coefficient in eq. ( |5.13| ). However, this coefficient is renormal- 
ization scheme dependent and requires special analysis (see below). 

In the scalar and pseudoscalar channels the general expression for the coefficient 



functions eq. ( p.45| ) takes the form 
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where 



m 



l/iTif is the quark mass renormalization constant (see eq. ( |2.25D ). The 
7^ matrix is defined within the dimensional regularization according to eq. ( |4.10| ). It is 
easy to see that in the calculations of Cq2, two matrices can be anticommuted over 
the fermion propagators and "annihilate" each other so that the results in both channels 
coincide. The calculational procedure is exactly the same as it was for the vector channel. 
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except for the need of mass renormalization. The results for the coefficient functions (7(^2 
and in the MS -scheme are as follows (Surguladze and Tkachov, 1986, 1990). 

In the (pseudo)scalar channel 



1 



Nplas 



TT 



In the scalar channel 



1 



TT 4 



a. 



TT 



C, 



/447 21 
V208 ~ 13^ 



m)+CA 



0{al) 



389 3 
T44 + 26«^) 



(5.16) 



5 25 

T TN 

39 36 



+ 0{al) 



and in the pseudoscalar channel 
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The result for coincides with the analogous one for the scalar channel. 

Let us turn to the renormalization group analysis of the above results, in this particu- 
lar case it is possible to use the following trick (Surguladze and Tkachov, 1990). Note first 
that the vacuum average of the rcnormalizcd operators and mqq and their coefficient 
functions depend on the renormalization parameter [i and therefore are not convenient 
for further analysis. However, as was shown by Collins, Duncan and Joglekar (1977) (see 
also Nielsen, 1977; Tarrach, 1982; Narison and Tarrach, 1983), the vacuum average of the 
trace of the energy-momentum tensor 
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is renormalization group invariant. On the other hand, in the MS type schemes the quark 
condensate < mfq^qj >o is renormalization group invariant to all orders of perturbation 
theory (see, e.g., Tarrach, 1982). One can introduce the renormalization group invariant 
quantity 
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so that the new coefficient functions defined from equation 



E < mfQfQf >o 
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CG4^,,as)n + CU^^,aAY.< mjqjqf >o (5.22) 
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should be the renormahzation group invariants. This is true since the l.h.s of eq. (|5.22| ) 
is directly connected to the observables (Shifman, Vainshtein and Zakharov, 1979) and 
consequently is invariant. From eqs. ( |5.21| ) and ( |5.22| ) we find the invariant coefficient 
functions corresponding to the invariant combinations of the gluon and quark condensates 
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Note that, in fact, there are terms of the type mjrrrij, or/and mj in the r.h.s. of eq. 
( ^.22| ). However, obviously these terms do not affect our equations for invariant coeffi- 
cient functions. The two-loop coefficient functions for ~ terms have been calculated 
by Chetyrkin, Gorishny and Spiridonov (1985). The contributions from such terms are 
negligible for phenomenological applications and will not be discussed here. 

Now one can use the renormahzation group invariance of the coefficient functions and 
write 

/ ,,2 \ 

12\ 



=C,(l,a,(Q^)). 



(5.25) 



Reevaluating the coefficient functions for the u, d, s hght quarks {N = 3) we obtain 
the following results in the MS scheme. 

In the vector channel 
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In the pseudoscalar channel 
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For all channels 
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(5.26) 
(5.27) 

(5.28) 
(5.29) 

(5.30) 
(5.31) 



Note again very large 0(q;^) corrections for the coefficient functions of quark condensates 
in the MS-schcmc. The running couphng is evaluated at the typical hadronic mass scale. 
Presently the 0{as) corrections have also been calculated for the dim = 6 operators 
(Lanin, Spiridonov, and Chetyrkin (1986)). We also mention the calculations in the case 
of heavy quark currents (see, e.g., Broadhurst et al, 1994 and references therein). 



VI R{s) in electron-positron annihilation to 0{al) 

In this section we present an outline of the evaluation of the corrections up to 0{al) to 
the total cross section in the process e+e~ hadrons (Fig. 5) in the limit of zero light 
quark masses and infinitely large top mass. We also mention the QCD evaluation of the 
hadronic decay rates of the Z boson and the relevant quark mass effects. 



FIG. 5. The process e+e^ hadrons. The shaded bulb inchides any interactions of 
quarks and gluons (or ghosts) allowed in QCD. The dots cover any relevant number of 

gluon and quark propagators. 

These calculations were first attempted by Gorishny, Kataev and Larin (1988). How- 
ever, it was shown that those results were incorrect. Indeed, about 4 years ago an indepen- 
dent calculation of the above quantity was completed (Surguladze and Samuel, 1991a,b). 
The result is much smaller and has the opposite sign compared with the old 1988 result. 
This finding was confirmed shortly after that by Gorishny, Kataev and Larin (1991). 

In the process shown in Fig. 5 an electron-positron pair annihilates producing either a 
photon or a Z-boson, which further produces quark-antiquark pairs (in QED) plus gluons 
(if strong interactions are "switched on"). Finally, quarks through hadronization form 
hadronic final states with probability equal to one (confinement hypothesis) and the total 
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cross-section is given by 



o-tot(e+e" hadrons) = — — 3 ^Q'j{l + 6qcu), (6.1) 

where s is the total centre-of-mass energy squared, Q/ is the electric charge of the partic- 
ipating at the given energy quark flavor /, factor 3 stands for the number of color degrees 
of freedom and 5qcd stands for the strong interaction contributions. The hadronic pro- 
duction in electron-positron annihilation is usually characterized in terms of the i?-ratio 
-the total hadronic cross section normalized by the muon pair production cross section 

G( \ ^tot(e+e- ^ hadrons) 

The above expressions are relevant at energies much less than the Z mass {^/s <^ Mz) 
corresponding to, for instance, PEP/PETRA energy range. At LEP the effects of the Z 
boson become important. The corresponding i?-ratio is defined as a ratio of the hadronic 
and electronic widths of the Z boson 

_ T{Z ^ hadrons) 

Note that the total hadronic width of the Z boson in the above equation is the sum of the 
vector and axial current induced decay rates. Strictly speaking, those rates get different 
strong interaction contributions. In the present section we calculate the QCD corrections 
in the vector channel - 5qco in the limit of massless light quarks and the infinitely large 
top mass. This quantity is, in fact, relevant for the axial part as well. To get the complete 
axial decay rate, additional contributions are necessary. For details see the original works: 
Kniehl and Kiihn (1990), Kniehl (1990), Chetyrkin and Kiihn (1990), Chetyrkin, Kiihn 
and Kwiatkowski (1992), Chetyrkin (1993a), Soper and Surguladze (1994), Surguladze 
(1994c), the review articles by Kniehl (1994b, 1995b), Soper and Surguladze (1995) and 
also section 4 of the present paper. 



A via renormalization constants 



The vacuum polarization function Il{Q^) defined in eq. ( ^.101 ) has a cut along the negative 
axis in the massless case. The ratio R{s) can be found taking the imaginary part of 



n(s + iO), according to eq. (|2.12|) . Alternatively, R{s) can also be found from eq. (|2.13|) , 



which in combination with eq. ( p.l6| ) gives to 0{ 



R(s) = fl, + ^fl, + c^Yr, + r^Vfe -^-^R.]. (6.4) 



71 Vvr/ \7r/V 3 

The origin of the large and negative scheme-scale independent term i?i7r^/3g/3 can be 
understood if one takes into account the presence of ~ log^/i^/s terms at 0{al) in the 
Il-f unction and 

-Imlog^(s + iO) = -31og^s -h TT^. 
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The leading QCD term Ri at O(a^) is due to the couphng renormahzation. Note, that the 
Ri in the above equation are the perturbative coefficients of the D{Q'^) function defined 
in eq. ( |2.13| ). For the discussion of the procedure of analytical continuation and the origin 
of additional ~ vr^ terms, see also Krasnikov and Pivovarov (1982), Pennington and Ross 
(1982), Radyushkin (1982), and Pivovarov (1992a). 

Substituting eq. (p.27|) with the renormalized strong coupling into eq. ( p.l2|) and 
taking into account the relations ( |2.29|) , (|2.3CI| ), (|2.32| ), ( p.33| ) we obtain 



R{s) = -l{Zi,-i + ^(2Z2,_i) + (3^3,-1 - m2,o + 2/?o^2,-ilog^ 



TT 



4Z4,_i - 2/3on3,o - An2,o + 2/32n2,i - ^^^2,-1 
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+ (6/3o^3,-i + 2AZ2,„i-2/5o%,o)log^ + 2/5o'Z2,_ilog2^] +0(«f)}. (6.5) 

Note that the appearance of perturbative coefficients of the renormahzation constant in 
the above equation is totally due to the relations ( p.33|) . In fact, Zn has only simple poles 
and hence has no imaginary part. The latter is the specific feature of the MS prescription. 
The expression ( |6.5| ) exhibits one of the main ideas of this calculation. Namely, in order to 
calculate the /-loop contribution to R, it suffices to calculate the /-loop counterterm Zu to 
the bare quantity 11^, and the / — I -loop approximation to 11^. In other words, the minimal 
information necessary to obtain the four-loop R{s) is contained in the divergent part of 
one-loop diagram, two-loop diagrams calculated up to ~ e terms, three-loop diagrams 
calculated up to the finite parts in the limit £ — >■ and only a leading ~ 1/e terms 
in the overall counterterms of the four- loop diagrams. In fact, as we demonstrate in the 
next subsection, using the infrared rearrangement procedure (Vladimirov, 1980; Chetyrkin 
and Tkachov, 1982) one can complete the entire calculation dealing effectively only with 
three-loop diagrams. We mention once again that, through the procedure of infrared 
rearangement, within the MS prescription, the problem of calculation of the counterterms 
to arbitrary / -loop diagrams with an arbitrary number of masses and external momenta 
can be reduced to the calculation of / — 1 -loop propagator type massless integrals up 
to finite terms. In our case / = 4. On the other hand, the recursive type algorithms 
for multiloop Feynman integrals (Chetyrkin and Tkachov, 1981; Tkachov, 1981, 1983) 
and their computer implementation (Surguladze and Tkachov, 1989; Surguladze, 1989b, c, 
1992; Gorishny, Larin, Surguladze and Tkachov, 1989) allow one to calculate propagator 
type Feynman diagrams to three-loop level. 



B Full calculational procedure with a typical four-loop diagram 

In this subsection we demonstrate the full calculational procedure for a typical four-loop 
diagram pictured in Fig. 6, which contributes to the photon renormahzation constant 
Zjj and hence to the i?-ratio. To simplify the description, in some cases we will avoid 
complicated equations, substituting for them their graphical representation. 
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FIG. 6. A typical four-loop nonplanar type diagram contributing to R{s) 



We need to evaluate the counterterm to the diagram pictured in Fig. 6. In other 
words, we should evaluate —ICR' for this diagram. A simple power counting shows that 
the given diagram diverges as 

G ~ lim Q^^-i^ 

and the superficial degree of divergence is 2. Using the fact that the counterterm has only 
a polynomial dependence on the external momenta Q within the MS prescription, one can 
remove such a dependence by differentiating the diagram twice with respect to Q and then 
set the external momentum to zero. At the next step, since there is no dependence on 
the external momentum, one can introduce a new fictitious external momentum flowing 
through one of the diagram fines. This line should be chosen in a way that simpfifies the 
topology of the diagram and avoids infrared divergences. The above procedure for the 
diagram in Fig. 6 is displayed in the following graphical equation 



Z D ]CR' 
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ICR'^^4{2-D) +2 I (6.6) 



= /C|4(2-D) +2 |, 
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where the dot and dashes on the hnes result from differentiating the corresponding fermion 
propagators 
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Boxes contain the corresponding three-loop propagator type subgraphs with subtracted 
divergences - complete /^-operation (Fig. 7). The dotted lines mean that this line is tem- 
porarily "torn" . After the evaluation of boxes, the parts of the torn line should be pasted 
and a trivial fourth loop integration should be done, taking into account the corresponding 
exponents of the propagators due to the three-, two- and one-loop "box" insertions. The 
above procedure gives a great simplification of the problem. Indeed, the evaluation of the 
four-loop counterterm is reduced to the evaluation of three-, two- and one-loop graphs. 



- A 



R 



B 



R 



FIG. 7. Complete i?-operation for the three-loop subgraphs 

The complete i?-operation of the three-loop diagram insertions corresponding to the 
ones at the r.h.s. of eq. ( |6.6|) is given in Fig. 7. Graphs in the brackets correspond to two- 
and three-loop counterterms. There is no one-loop divergent subgraph in this particular 
diagram. Thus, 

where Gi is any divergent subgraph of the given diagram. It is easy to recognize that the 
two-loop subgraph in Fig. 7 does not have sub divergences (only an overall one) and the 
corresponding counterterm is simply the pole part of this subgraph 
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ICR'l = /C 
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Analogously, because of the topology, the three-loop counterterm does not have a 
subdivergence and the corresponding counterterm is the pole part of this diagram 



]CR'\ = /C 



If, in general, a diagram contains divergent subgraphs, then the recursive formula 
( p.35| ) should be used. 

As a result of the above manipulations, we managed to reduce the problem of calcu- 
lation of the counterterm to the four-loop diagram pictured in Fig. 6 to the calculation 
of several three-, two- and one-loop diagrams shown in Fig. 7. Note, however that the 
"dots" and "dashes" on the diagram lines make their evaluation significantly more dif- 
ficult. The computer programs for analytical programming systems capable of handling 
such calculations are the SCHOONSCHIP program MINCER (Gorishny, Larin, Surguladze 
and Tkachov, 1989; Surguladze, 1989b, c) and the FORM program HEPLoops (Surguladze, 
1992). The latter is especially well-suited for large scale calculations and is much more 
efficient than the MINCER program. 

It is important to stress that, in fact, it is sufficient to evaluate only the ICR' for the 
relevant three-loop subgraphs. In other words, it is not necessary to calculate separately 
three-loop counterterms similar to the graph in the last brackets for the box A in Fig. 7. 
Indeed, a more detailed analysis gives 

R[G] = R'[G] - (1 - D/2)}C {jrr^R'[G]j , (6-7) 

where G is the corresponding three-loop subgraph. The above relation allows simple 
computer implementation and facilitates calculations considerably. 

The complete i?-operation for each four-loop diagram generally has the form 

(2 \ 4£ 3 / 2 \ (4—0^ 

where fi{e) is the result of the calculation of the corresponding Feynman graphs including 
the last trivial loop integration and q are the /-loop counterterms, polynomials in l/e. 
As we already mentioned, in the MS type renormalization scheme, the counterterm for 
a particular diagram is a polynomial in dimensional parameters (see, e.g., the textbook 
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by Collins, 1984 and references therein). Thus, the terms of the type (l/s)" ln"^(yU^/(5^), 
which appear due to the expansion of the factors {^'^/Q'^Y^ into the Laurent series in e, 
must be canceled in the final answer for the particular diagram. This can be used to test 
the calculations at the graph-by-graph level. Recall, that we calculate the counterterm 
Z4 j to the four-loop diagram. 

Finally, for the contribution to the of the diagram pictured in Fig. 6 we obtain 
the following result 

(g)' 7V.C-,(C.-C.)(C.-C./2) ^ [4-26J, + f i-40C(3)l;. 

The CPU time for the above diagram on a 0.8 MFlop IBM compatible mainframe was over 
6 hours. The extended version of the program MINCER for the system SCHOONSCHIP was 

used. Note that the above result, as well as the total result for the photon rcnormalization 
constant does not depend on any modification of the minimal subtraction prescription. 



C Four- loop results 

In this subsection, we present results and some of the details of the 0{al) QCD evaluation 
of the ratio R{s) in electron-positron annihilation (Surguladze and Samuel, 1991a,b). 

The total number of topologically distinct Feynman diagrams contributing to Zi ^ is 
1, to Z2^i is 2, to 4 is 17 and to Z4 j is 98. However, after application of the infrared 
rearrangement procedure which, as discussed above, involves differentiation twice with 
respect to the external momentum of the diagram, the number of four-loop graphs which 
need to be calculated increases to approximately 250. Furthermore, there are one-, two- 
and three-loop diagrams, approximately 600, which need to be calculated to subtract 
sub divergences (evaluate R') for all four-loop diagrams. 

All analytical calculations of the four-loop diagrams have been done by using the 
program, which is an extended version (Surguladze, 1989c) of the program MINCER (Gor- 
ishny, Larin, Surguladze and Tkachov, 1989; Surguladze, 1989b). This version includes 
new subprograms for 4th loop integration and for ultraviolet rcnormalization. Evalua- 
tion of one- and two-loop counterterms has been done by using the program LOOPS 
(Surguladze and Tkachov, 1989a). The above programs are written on the algebraic pro- 
gramming systems SCHOONSCHIP (Veltman, 1967; Strubbe, 1974) and REDUCE (Hearn, 
1973) respectively. The full calculation took over 700 hours of CPU time on three IBM 
compatible 0.8 MFlop EC-1037 mainframes with the SCHOONSCHIP system. We have 
also recalculated some of the difficult four-loop diagrams with HEPLoops - a new program 
for analytical multiloop calculations (Surguladze, 1992). The status of these and some 
other programs has been reviewed recently in Surguladze (1994d). 

In the diagram calculations the Feynman gauge is used. The momentum integrations 
are performed within the MS modification (Bardeen, Buras, Duke and Muta, 1978) of 
the minimal subtraction prescription ('t Hooft, 1973), which amounts to formally setting 
7 = C(2) = log47r = 0. A discussion of the scheme dependence of the results is given at 
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the end of this section and in section 9. The full graph-by-graph results will be published 
elsewhere. 

The analytical result for the four-loop photon renormalization constant reads 



Zph = 1 + -^Zu — 

£2^^ 9 ^ 162 ^ 9 ^81 81 / 

It should be stressed that the Riemann (^-functions ({4) and C(5), which appear at the 
individual graph level cancel in the above expression. Moreover, as we have observed, 
the C(4) has disappeared within each gauge invariant set of diagrams. Note that C(3) 
disappears for QED {Cp = 1, Ca = 0, T = 1) except the last term, which comes from the 
"light-by-light" type diagrams (Fig. 8). The diagrams pictured in Fig. 8 are some of the 
most complicated ones and the computation of each of them requires over 80h of CPU time. 
Note, however, that the second and fourth diagrams in Fig. 8 differ correspondingly from 
the first and third ones only by the SU(N) group weights. So, in fact, only two of them 
have been calculated. The result ( |6.8| ) does not depend on the particular modification of 
the minimal subtraction prescription. 



FIG. 8. "Light-by-hght" type diagrams 
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In order to evaluate R{s) to 0{a^), besides the four-loop we calculate the un- 
renormalized hadronic vacuum polarization function Il^{Q'^) to the three- loop level. We 
get the following analytical result in the MS scheme. 
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(6.9) 



The above result depends on the particular modifications of the minimal subtraction 
prescription, unlike the result for the renormalization constant ( |6.8| ). 

Substituting the expressions for the relevant Zjj- and Iljj-, extracted by comparing 
eqs. ( ^TsD and (|6.9|) to eqs. (|2.26| ) and ( p.27|) , into eq. (|6.5|) and recalling the values for (3q 
and Pi from eq. ( p.8|) we get the following MS analytical result for R{s) at the four-loop 
level 
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(6.10) 
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The logarithmic contributions are absorbed in the running couphng by taking /i^ = s. 
Those contributions will be presented explicitly in section 9. Note that C(5) appears in 
the final result due to the contributions from Hsfi. The last term ~ iJlfQfY comes from 
the so called "light-by-light" type diagrams (Fig. 8). For standard QCD with the SUc(3) 
gauge group we obtain 
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(6.11) 



Finally, taking into account the values for the relevant Ricmann -functions, ^(2) = 7r^/6, 
C(3) = 1.2020569... and C(5) = 1.0369278... we obtain the numerical form 
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TT 
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TT 
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(6.12) 



Note that only 19 four-loop diagrams contribute to the term ~ N and 2 four-loop diagrams 
contribute to the term ~ N^. The most comphcated diagrams are pictured in Fig. 9. The 
CPU time for each of them was over lOOh and the intermediate expression had as many 
as ~ 10^ - 10^ terms. 



FIG. 9. Some of the most complicated diagrams 
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It is known, that the perturbative coefficients for R{s) are scheme dependent. The 
above resuh was obtained in the modified minimal subtraction, the so-called MS scheme 
introduced by Bardeen, Buras, Duke and Muta (1978). While the scheme-scale dependence 
problem will be discussed in section 9, here we present the results for a couple of other 
versions of the minimal subtraction scheme. First, we consider the so called G scheme 
(Chetyrkin and Tkachov, 1979, 1981; Chetyrkin, Kataev and Tkachov, 1980), which is 
convenient for practical multiloop calculations. The G scheme is defined in such a way 
that the trivial one-loop integral in this scheme is 
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The result for R{s) in this scheme is 
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(6.13) 



The parametrization of the running coupling in the above equation has the same form as 
in eq. ( |2.17D . However, the parameter A has to be changed to some other parameter A^. 
Finally, in the original MS scheme ('t Hooft, 1973) we get 
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(6.14) 



As one can see, starting from O(a^) the results heavily depend on the choice of 
the particular modifications of the minimal subtraction scheme. This dependence, called 
renormalization group ambiguity of perturbative results is an important problem and 
deserves special consideration. We will return to this issue in section 9. 

Concluding this section, we mention once again that the results of the above described 
calculation of the four-loop correction to the -R(s) have been published in Surguladze and 
Samuel (1991a,b) and independently Q in Gorishny, Kataev and Larin (1991) and hence, 

•^See however the discussion in the last three paragraphs of section 3 in the review article by Surguladze 
(1994d). 
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most likely, the above results are reliable. Interesting relations between the radiative 
corrections for different observables, found by Brodsky and Lu (1994, 1995) serve, in 
particular, as another confirmation of our results. 



VII r(r- ^Vr + hadrons) to 0(af) 

The other important inclusive process for phenomenology and testing the Standard Model 
is the hadronic decay of the r lepton (Fig. 10). For a recent review see, for instance, Pich 
(1994a). For earlier references see AltareUi (1992), Marciano (1992), and Pich (1991). 



FIG. 10. Hadronic decay of the r- lepton 

In this section, using our result of four-loop calculation of the (Ttot(e''"e~ — > hadrons) 
(Surguladze and Samuel, 1991a,b), wc evaluate the hadronic decay rate of the r lepton to 
O(a^) in pcrturbativc QCD (Pich, 1990; Gorishny, Kataev and Larin, 1991; Samuel and 
Surguladze, 1991; see also Braaten, Narison and Pich, 1991; Pich, 1992a,b; Diberder and 
Pich, 1992a,b; Pivovarov, 1992b). We also comment on the status of the nonperturbative 
corrections to this quantity. 

We follow the method first suggested by Tsai (1971), Shankar (1977), and Lam and 
Yan (1977) for theoretical evaluation of heavy lepton decay rates. This method has been 
further developed for the r lepton including the higher order pcrturbative corrections 
and involving the operator product expansion technique (Wilson, 1969) to analyze the 
nonperturbative contributions (Schilcher and Tran, 1984; Braaten, 1988; Narison and 
Pich, 1988). As was shown in the above works, combining the operator product expansion 
technique and analyticity properties of the correlation function of quark currents, the ratio 

V{t- + hadrons) 

is calculable in pcrturbative QCD. Strictly speaking, besides the QCD pcrturbative parts 
the nonperturbative and weak contributions should be included to estimate Rt- There 
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are instanton contributions as well. However, it was shown recently by Nason and Por- 
rati (1993) (see also Kartvelishvili and Margvelashvili, 1995) that these contributions are 
completely negligible due to the chiral suppression factor uiuTniims/ M"^. The can be 
written as the following sum 

Rr = Rr' + + (7.2) 



A Perturbative QCD contributions 



The quantity can be expressed as the following integral over the invariant mass of 
the hadronic decay products of the r lepton (Lam and Yan, 1977; Braaten, 1988) 



^pert 



3 jM-^ ds / s y 



1 + 2-^)lmn^{s + iO) + Imn^(s + iO) 



(7.3) 



where is the mass of the r lepton. The functions 11^ and 11^ are the transverse 
and longitudinal parts of the correlation function of weak currents of quarks coupled 
to W boson. In fact, 11^'^ are the appropriate combinations of vector and axial parts 
corresponding to the vector and axial currents of u, d, s light quarks (for details see, e.g., 
Pich, 1994a). The expression for in the form of ( [7.3|) is not quite useful. The problem 
is that the correlation functions involved can not be calculated at low energies because of 
the large nonperturbative effects that invalidate perturbative approach. However, simple 
analyticity properties of the correlation functions allow us to evaluate the integral in (773). 
Indeed, the function H is analytic in the complex s plane everywhere except the positive 
real axis. According to the Cauchy integral theorem, an integral over s along the closed 
contour Ci + C2 (Fig. 11) of the product of H(s) with any nonsingular function /(s) is 
zero. 



FIG. 11. Integration contour 
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On the other hand, the imaginary part of the correlation function is proportional to 
its discontinuity across the positive real axis. So, the following relation holds 



Mi 



ds /(s)Imn(s) 



i/ dsfisMs) = -^[ dsfisMs) 



(7.4) 



where the C2 is the circle of radius | s |= (Fig. 11). The benefit of the above relation 
is that in the r.h.s. one needs to calculate the correlation function for | s | at M^. 
Hopefully, is large enough to use the operator product expansion in powers of 1/M^ 
and the q;s(M^) is small enough to use perturbative expansion in a^. Then the perturbative 
method can, in principle, be used to calculate the leading term in the operator product 
expansion and the higher twist terms can be estimated semi-phenomenologically. 

Using eq. (|7.4| ), the perturbative part of the ratio can be expressed by an integral 
over the invariant mass s of the final state hadrons along the contour C2 in the complex 
s-plane (Fig. 11). In the chiral limit, = = = 0, the currents are conserved and 
the longitudinal part of the n(s) is absent. In the axial channel n'^(s) = 0{m'j/s) (see 
section 4). For the R^"^^^ we get 



^pert 



3i r ds 
8^Jc2 



1--V 
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'M2 



(7.5) 



Note that the factor (1 - s/M^f suppresses the contribution from the region near the 
positive real axis where the n(s) has a branch cut (Braaten, 1988). To simplify the 
description, we use the chiral limit which is a perfect approximation for Rr. On the other 
hand, the mass corrections can be included with the calculation very similar to that in 
section 4. The actual calculations show (Chetyrkin and Kwiatkowski, 1993; see also recent 
analyses in Pich, 1994a) that the effects of quark mass corrections on i?^ are well below 
1% and can be neglected. Note also that, in the massless quark limit the contributions 
from vector and axial channels to 11 coincide at any given order of perturbation theory and 
evidently the results are flavor independent. So, in this case, for evaluation of n"^(s) in 



eq. ( |7.5| ) we use our earlier results for the electromagnetic two-point correlation function 
that contributes to R(s) in electron-positron annihilation (section 6). 

The function Il^(s) can be related to the D(s) function defined in section 2 as follows 



3 d 
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U'^is) 



f=d,s I K/ 



-D( 



(7.6) 



where Vud and V^s are the Kobayashi-Maskawa matrix elements. | Vud P + I Vus P= 
0.998 ± 0.002 (see, e.g., Pich, 1994b). The factor in the r.h.s of eq. (|^) is due to the 
replacement of the electromagnetic currents by charged weak currents in the correlation 
function. Note also that evidently the "light-by-light" type graphs (Fig. 8) do not con- 
tribute to the decay width of the r lepton. Thus, the term ~ {J^fQfY drops out in the 
D function. The perturbative coefficients of D{s) have been given in the previous section 
up to the four- loop level in the vector channel (see eqs. (|2.16|) and (|6.11| ) ). 
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Performing the contour integration in eq. ( [7.51) using the relations ( |7.6| ) and ( |2.16|) , 
and replacing as{s) by ttsl-^r) using the evolution equation ( |3.19| ), we obtain in terms of 
perturbative coefficients of R{s) 
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where, as we have already mentioned, the term ~ (Z)/ QfY should be omitted in R^. 

Substituting the relevant expressions for Ri and Pi from the previous sections, we 
obtain the 0{a^) analytical result in the MS scheme 
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Within the standard QCD with the SUc(3) gauge group we obtain 
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and a numerical form reads 
3(0.998 ±0.002) 



1 + + 5.2023 



26.366 
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TT 



TT 



(7.10) 



B On the Nonperturbative and Electroweak contributions 

The nonperturbative contributions to Rr can be expressed as a power series of corrections 
in 1/M2 

ononpert ^2 { (^^ (M.) , C^. < O2. >0 

where the m/ are u, d, s running quark masses, < >o are the so-called vacuum conden- 
sates, which can be obtained phenomenologically and the Ci are their coefficient functions 
describing short distance effects. Note that, in eq. (|7.11| ) we formally include part of 



the pure perturbative corrections (the first term) which is due to the nonvanishing m, d, s 
quark masses. These corrections for the u and d quarks are completely negligible. The 
contribution coming from the s quark is suppressed by sin^ 9c and is also below 1% (Pich, 
1990). Presently, the only way to estimate the strong interaction effects in the condensate 
contributions is by perturbation theory. The coefficient functions are asymptotic per- 
turbative series in terms of as- In order to estimate the nonperturbative contributions, one 
needs to sum up the power series of the QCD perturbative series. In the previous section 
we have described the calculation of the high-order perturbative QCD contributions to the 
coefficient functions of the dimension 4 power corrections (gluon, < a^G^ >o and quark, 
< TTLfqjqf >o condensates). It was shown (Loladze, Surguladze and Tkachov, 1985; Sur- 
guladze and Tkachov, 1989b, 1990) that the high-order perturbative corrections to some 
of the coefficient functions are too large. For instance, for the coefficient function of the 
condensate < m^ss >o in the vector channel (see eq. ( |5.27| )) Aefr ~ 30Aj^. This indicates 
that the renormalization group invariant criteria to the perturbative calculability of the 
QCD contributions to the coefficient function is not fulfilled. The coefficient functions of 
the dimension 6 condensates are calculated up to 0{as) (Lanin, Spiridonov and Chetyrkin, 
1986) and to analyze the corresponding series one needs at least the next to leading cor- 
rection. The above uncertainty in coefficient functions C2i allows one to estimate the 
condensate contributions probably not better than their order of magnitude. There is 
another source of theoretical uncertainties in the evaluation of condensate contributions 
of dimension 6 and higher, where the operator basis of expansion includes a large number 
of operators. Presently, there are no precise methods to estimate their matrix elements. 
For the matrix elements of four quark operators (dimension 6) the vacuum saturation 
approximation (Shifman, Vainshtein and Zakharov, 1979) is used to express them as the 
square of the two-quark matrix elements. However, the vacuum saturation approximation 
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is not expected to be precise enough in order to use it in the analyses of the tiny non- 
perturbative contributions (see, e.g., analysis by Altarelli, 1992; see also a brief discussion 
in Surguladze and Samuel, 1992b ). Indeed, as it was found by Braaten (1988) and Pich 
(1990, 1992a,b, 1994a), the nonperturbative corrections are below the 1% level with large 
theoretical error. The contributions of dim=4 condensates start at 0{a1) and thus are 
suppressed by two powers of as- The dim=6 and dim=8 corrections are suppressed by the 
inverse powers of Mr (M^ and respectively) and are small. On the other hand, the 
corrections in vector and axial channels have opposite signs and they largely cancel each 
other, so the total relative error is even larger. In the works by Pumplin (1989, 1990) it 
was shown that the uncertainty due to threshold effects makes a significant contribution 
in the theoretical error for i?^. In the works by Altarelli (1992) and Altarelli, Nason and 
Ridolfi (1994) an ambiguity ~ A^/M^ is discussed. Earlier, Zakharov (1992) has argued 
that such dim=2 terms in eq. ( [7.11| ) can be generated by ultraviolet renormalons. For 
an alternative point of view on the effects of possible dim=2 terms, see Narison (1994). 
However, this issue is still a subject of intensive discussions and likely is far from being 
settled. 

Summarizing, we note that the above mentioned major sources of theoretical un- 
certainties in the evaluation of small power corrections makes certain restriction on the 
precision theoretical prediction of Rr and consequently on as(M^). Fortunately, the non- 
perturbative corrections are suppressed and the hadronic decay of the r still remains as a 
good source to extract the low energy a^. 

Finally, we note that the electroweak contributions R^'^^^ were calculated by Marciano 
and Sirlin (1988), and Braaten and Li (1990). Those corrections contain logarithms of 
Mt/Mz and are not negligible. The leading order electroweak corrections give roughly 
+2% contributions to Rr (see, e.g., Pich, 1994a). 



VIII Four-loop QED Renormalization Group Functions 

In this section we outline the calculation of the standard QED renormalization group func- 
tions at the four-loop level in the minimal and momentum subtraction schemes. These 
quantities can be obtained as an intermediate result of the calculations of R{s), described 
in the previous sections, by replacing the SUc(3) gauge group invariants for the corre- 
sponding diagrams in a proper way. The results of two independent calculations of the 
four-loop QED /3-function by Gorishny, Kataev and Larin (1990), and by Surguladze 
(1990) have been reported in the joint publications by Gorishny, Kataev, Larin and Sur- 
guladze (1991a,c). 
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A General formulae 

The Lagrangian density of standard QED is 
Lqed = 

--F^.F^^ + tY^i^^rD,i^j - E^i^iV^. - T^d.A'^d.A'^, (8.1) 

where F^^, = d^A^ — d^A^ and = — ieA^. is the gauge parameter, nij are the 
fermion masses, ip and are the fermion and photon fields and e is the electric charge. 

Renormalization constants are defined by the relations 



= ^-yZphA^ (8.2) 



For the fermion-fermion-photon vertex renormalization one has 

/i-2^Z.erte?7M^^V^ = ^i-^'Jz^Z^e^^^A''^p. (8.3) 



According to Ward identity in QED (Ward, 1950) Z^crt = Zp, which implies from eq. (|8l3| ) 
the identity 

ZaZp^ = 1. (8.4) 

From eqs. (|8.2| ) and (|8.4|) we get 

ae = f^^'Z^^a. (8.5) 

The gauge invariance of the QED lagrangian implies the absence of the conterterm for the 
gauge fixing term in (|8.1|) and, thus, Zq = Zp^. 



Using the relation ( p.5|) and the renormalization group invariance of "bare" coupling 
fi^daj^/dfi"^ = 0, taking into account that Z^^ depends on /i only via a and also the 
standard definition of the QED MS /3-function 



(8.6) 

qb fixed 



we obtain a convenient expression for the further evaluation of the /3 function 

47r £^0 1 - a^logZph 
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B Four-loop results 



The photon field renormahzation constant can be found from the QED relation, anal- 



ogous to eq. (|6.8|) , where only 58 QED four-loop diagrams contribute to n(/i /Q , a). The 
prescription for the evaluation of the diagram contributions to the IIb is analogous to the 
one described in section 2. The total CPU time on the three IBM compatible mainframes 
was approximately 400 hours. Setting Cp = 1, Ca = 0, T = 1 and = a in eq. (|6.8|) , 
we obtain the four- loop photon renormahzation constant in QED, corresponding to the 
minimal subtraction prescription 
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Substituting the expression for Zph into eq. (p.7|), we obtain the following result for the 
four-loop QED /3-function in the MS type schemes 
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It is useful for further applications to present the result for the Johnson- Willey-Baker 
Fi function (Johnson, Willey and Baker, 1967; Baker and Johnson, 1971; Johnson and 
Baker, 1973). This function can be obtained from the result for P^^^ by subtracting 
the contributions of the diagrams with fermion loop insertions into the photon lines and 
reducing the power in a /An by one. We obtain 

Note that all coefficients up to four-loop level are rational numbers. The results for most 
of the individual graphs do contain transcendental C(3), C(4) and C(5)- The C(4) and C(5) 
cancel within each gauge-invariant set of diagrams. The three-loop results agree with the 
ones obtained by de Rafael and Rosner (1974). It is possible to recalculate the MS QED 
f3 function in the form of the Gell-Man-Low \I'(q;) function - the QED (3 function in the 
MOM scheme. See details in Gorishny, Kataev, Larin and Surguladze (1991a) (see also 
Adler, 1972; de Rafael and Rosner, 1974). We obtain the Gell-Mann-Low \1/ function at 
the four-loop level 
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The 0(q;^) result agrees with the one obtained by Baker and Johnson (1969) and Acharya 

and Nigam (1978, 1985). 

Recently, Broadhurst, Kataev and Tarasov (1993) have carried out an additional cal- 
culation necessary to convert the four-loop MS QED {3 function to the four-loop QED 
on-shell (3 function, usually called the Callan-Symanzik function fi"^^ (Callan, 1970; 
Symanzik, 1970, 1971). This function is defined as follows 



(8.12) 

ofB fixed 



where rrie is the electron pole mass. The subtraction prescription in this case requires all 
subtractions to be on-shell. The three-loop Pq^d calculated long ago by de Rafael 
and Rosner (1974). The four-loop result has the following form (Broadhurst, Kataev and 
Tarasov, 1993) 



i-(l->).?a.)..-ic(3)).-(gl->-». 



a' 
(8.13) 



IX Renormalization Group Ambiguity of Perturbative QCD 
Predictions 

In the previous sections we have demonstrated the calculation of some of the important 
observables within the framework of perturbative QCD. This involves calculation of a 
large number of Feynman diagrams and requires a a very large amount of computer and 
human resources. For example, to 0{ag) we have calculated 98 (effectively 250) four-loop 
Feynman diagrams. The next order requires calculation of approximately 600-700 five- 
loop diagrams. Calculations of such a scale are extremely difficult. On the other hand, 
perturbative QCD scries arc asymptotic ones and the question of how many orders need to 
be calculated, can be answered only from estimates of remainders (see, e.g., the textbook 
by CoUins, 1984). Moreover, perturbative coefficients beyond the two- loop level, as well 
as the expansion parameter, are scheme-scale dependent. The scheme-scale ambiguity - a 
fundamental property of the renormalization group calculations in QCD, does not allow 
one to obtain reliable estimates from the first few calculated terms without involving 
additional criteria. 
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In this section we discuss the extraction of rehable estimates for observable quantities 
within perturbation theory. The problem of scheme-scale dependence of perturbative QCD 
predictions will be considered first within the MS prescription and then we outline a scheme 
invariant approach along the lines of Stevenson (1981a,b). We apply the three known 
approaches for resolving the scheme-scale ambiguity. As a result, we fix the scheme-scale 
parameter, within the framework of MS prescription, for which all of the criteria tested are 
satisfied for the quantity R{s) at the four-loop level (Surguladze and Samuel, 1993). On 
the other hand, we estimate the theoretical error by using the scheme-scale dependence 
as a measure of the theoretical uncertainty (Surguladze and Samuel, 1993; Surguladze, 
1994b). We also mention the recent discovery of commensurate scale relations by Brodsky 
and Lu (1994, 1995). These relations allow one to connect several physical observables, 
providing important tests of QCD without scheme-scale ambiguity. 



A Perturbative QCD series: How many loops should be evaluated? 

The R-ratio in electron-positron annihilation is given within perturbation theory in the 
following form 

\ yU"^ TT /i"^ TT fi"^ IT ) 

Our further discussion is quite general and can be applied to other observables like or 
Higgs decay rates. We consider high enough energies, where R is a function of a single 
variable - the center-of-mass energy squared. Our aim is to evaluate pure QCD effects in R, 
which start with the term 0{a^^ within the minimal subtraction prescription ('t Hooft, 
1973). We should stress here that the calculational methods allowing one to evaluate 
perturbative corrections up to the four-loop order (up to the five-loop in some cases) is 
essentially based on some of the unique features of the MS prescription and our choice 
seems to be well justified. There is an ambiguity in the choice of renormalization scale 
parameter [i. Usually we set ^ = s and absorb the large logarithms in the definition of the 
running coupling. On the other hand, the choice = (x = e*) for all x gives equivalent 
expansions. Evidently, the sum of "all" terms in eq. ( |9.1| ) does not depend on the choice 
of /i. However, in practice, we deal with truncated series, where the sum has a nontrivial 
dependence on the choice of renormalization parameter. Here we keep the "natural" choice 
jj^ = s and the ambiguity is transferred to the prescription / d'^p — s> / c/^~^^p(/i^e*^*"^'^''^-*'')^. 
By changing t one gets different MS type schemes. One can always reexpand ( p.l[ ) in a 
new scheme (with a new A in ( |2.17D ) and so redistribute the values of (i > 1). All these 
schemes are equivalent. On the other hand, a new scheme may be "better", but one can 
conclude this only based on the knowledge of remainders. The problem of scheme-scale 
ambiguity which, in fact, is a problem of remainders can be formulated as follows. How 
does one choose ("optimize") the scheme (or A) in order to make the remainder minimal 
in the series of the type ( ^g. i|j for the given range of energy and what is the numerical 



69 



uncertainty of the approximation (9A)? Here one should also distinguish the following 
two questions. First, what is the best accuracy to which the given quantity is calculable 
via perturbation theory? Second, what is the accuracy of the given approximation? A 
few notes are in order. It is known that perturbative QCD series are asymptotic ones. No 
reliable estimates of the remainders are known at present. However, it is known from the 
theory of asymptotic series (see, e.g.. Dingle, 1973) that 

M 

\Y.ria\s)-R{s)\=Rx^AR^i^, when AT^A^opt. (9.2) 

i=l 

This means that, the remainder Rj^ goes to its minimal value AR^^^^ when the number of 
orders goes to its optimal value A^opt. Inclusion of the next to A/'opt orders will lead away 
from the correct value. It is known (see, e.g.. Dingle, 1973) that for a sign-alternating 
asymptotic series the remainder can be estimated by the first neglected term (or by the 
last included term). However, it is still unknown if the QCD perturbative series has this 
character. We assume as a hypothesis that within QCD one can estimate the remainder 
by the first neglected or last included term. Now, the minimal possible error, which de- 
fines the best accuracy of the perturbation theory for the given quantity has an order of 
^-Rmin ~ fx^ia-'^~^^{s), M — A/'opt. Notc that, both the number A/'opt and the value of 
the Ai?inin depend on the range of energy for the given process. We once again emphasize 
that the remainder depends on the choice of particular scheme and scale parameters and 
its estimate makes sense only for the "optimized" renormalization scheme which is unique 
for the given physical observable. In fact, it was argued (Stevenson, 1984, 1994) that, 
the "optimized" series can still converge even when the series in any fixed renormaliza- 
tion scheme is factorially divergent, if the "optimized" couplant shrinks in higher orders 
(see also Buckley, Duncan and Jones, 1993). However, whether this applies to QCD is 
unknown. 



B R{s) within the one parametric family of the MS type schemes and scale ambiguity 
problem 

Using the results of our four-loop calculations, we obtain the analytical result for R{s) with 
perturbative coefficients explicitly depending on the scheme-scale parameter (Surguladze 
and Samuel, 1993) 

Ris,t) = R^ + ^iMi?^ + (^A^]\R^ + f3^R^t) 

+ (^^)'[^3 - y/^o'^i + + PiRi)t + P'oRA- 

(9.3) 



70 



Recalling the values of the MS perturbative coefficients Ri from eqs. ( |6.4| ) and ( |6.10 ) and 
the (3i coefficients from eq. (|2.8|), we obtain numerically 



R(s, t) = 3 5] gjjl + + (^^i^)' [(1.9857 + 2.75t) - iV(0.1153 + 0.1667t)] 



+ ( ^^^^^^[(-6.6369 + 17.2964t+ 7.5625t2) 



TC 



-A^(1.2001 + 2.0877t + 0.9167t^) 
+Ar2(-0.0052 + 0.0384t + 0.0278^^)] 

-(E^/l f^^)'l.2395 + 0(a^), (9.4) 



/ 



where as{s,t) can be parametrized in the form of ( |2.17| ) with fi = s and A ^ = 
e~*/^Aj^. Obviously, t = corresponds to the MS scheme ( eq. ( |6.12D ). t = ln47r — 7 
will transform the result to the original MS scheme ('t Hooft, 1973). ( eq. ( |6.14D ). 
t=-2 corresponds to the G scheme introduced by Chetyrkin and Tkachov (1979, 1981) ( 
eq. (|6.13|) ). Note that because of a one-parametric nature of the MS prescription, the t- 
dependent terms in eq. (|9.4|) would represent also the scale dependence of the perturbative 
coefficients within the MS if one changes t log /i^/s and takes as{s,t) with s replaced 
by yU^ and t = 0. 

Several approaches were suggested to deal with the scheme- scale-remainder problem. 
Among them we consider the following ones. Fastest Apparent Convergence (FAC) (Grun- 
berg, 1980, 1982, 1984), where the next to leading perturbative correction is absorbed in 
the definition of the "effective" running coupling and the scheme-scale parameter is fixed 
accordingly. Principle of Minimal Sensitivity (PMS) of the approximant to the varia- 
tion of nonphysical parameters (Stevenson, 1981a,b, 1982, 1984; see also Mattingly and 
Stevenson, 1992, 1994). Brodsky- Lepage- Mackenzie (BLM) approach (Brodsky, Lepage 
and Mackenzie, 1983), which suggests one fix the scale by the size of the quark vacuum 
polarization effects resulting in the independence of the next to leading order perturbative 
correction of the number of quark flavors A^. For discussions of the above scheme-scale set- 
ting methods see Celmaster and Stevenson (1983), Brodsky and Lu (1992), and Stevenson 
(1992). The optimization of perturbation theory has previously been studied by Kramer 
and Lampe (1988), and Bethke (1989) for jet cross sections in electron positron annihila- 
tion. The optimized perturbation theory is tested for different physical quantities in QED 
and QCD by Field (1993). The scale ambiguity problem has been considered by Lu and 
de Melo (1991) for the (p^ model. The scheme-scale ambiguity problem for the quantities 
R{s) and Rr has been discussed by Maxwell and Nicholls (1990), Chyla, Kataev and Larin 
(1991), and Grunberg and Kataev (1992). Further study of the PMS method has been 
done in Raczka (1995). 

We apply the above methods to eq. ( |9.1|) and we find a scale which gives good results 
for all criteria considered (Surguladze and Samuel, 1993). We start by noting that, in gen- 
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eral, the renormalizatin scheme-scale dependence of perturbative resuhs are parametrized 
by the scale parameter, say, fi and the renormalization prescription dependent coefficients 
of P function (Stevenson, 1981a,b). We should stress however, that the /3 function is 
independent of any modification of the MS type prescriptions, but starting from /?2, the 
coefficients of (3 function do depend on the particular choice of subtraction prescription 
other than MS. In order to better visualize our discussion, we consider ffist the optimiza- 
tion procedures within the MS prescription. In other words, we fix the scheme dependent 
perturbative coefficients of f3 function to their MS values and consider only the scale 
variation. 



In Fig. 12 we have plotted r3(t) for different (see eqs. ( |9.1|) and ( |9.4| )). As one can 
see, within the region t ~ (—1.5, —0.5) has a very weak dependence on the number of 
flavors as well as on the parameter t. 



FIG. 12. r3(t) for different A^ 



Corresponding to the three-loop coefficient r2{t), straight lines intersect in one point for t ~ 
—0.7, which is obvious from eq. (|9.4| ). This value corresponds to the BLM result (Brodsky, 
Lepage and Mackenzie, 1983) /i^ = /x|jge°'^^°, and at this scale the flavor dependence is 
absorbed into the deflnition of the coupling. 

In Fig. 13 we have plotted the dependence of the partial sums 

n 

Rn{t) = E rUt){as/7cr, n = 1, 2, 3 

m=l 



on the parameter t. Here the parametrisation (|2.17|) was used, log s/A^ = 9 and A^ = 5. 



The general picture does not change for other reasonable values of log and A^. One can 



72 



see that PMS (Stevenson, 1981) works perfectly for a wide range of the logarithmic scale 
parameter t ~ (—1, +3) for the four-loop approximant and t ~ (—2, 0) for the three-loop 
approximant. A similar analysis at the three-loop level was done by Radyushkin (1983). 
According to the above analysis we found that the BLM scale t = —0.710 is good at 
the four-loop level as well (Fig. 12) and this value is within minimal sensitivity region 
(Fig. 13). 



FIG. 13. The approximants i?„ vs the scale parameter t 



Moreover, we found that if the t-parameter is chosen in the following jmalytical form 
t = 4:({3) — ll/2 + 0{e), which is equivalent to the definition of a new, say, MS modification 
of the MS scheme 

A~ = exp[-2C(3) + 11/4 + 0(e)]A^, (9.5) 

then the dependence and the C(3) terms cancel exactly at the 3-loop level. As a result, 
r2 = 1/12. Within this scheme the four-loop correction is almost independent of the 
number of flavors. The full result for the R-ratio for the arbitrary number of flavors can 
be written in the following simple form 



TT 12 V TT 



— ) (16.2 ±0.5) 

TT 



/ 



1.2 + 0(a^ 



(9.6) 

where the small uncertainty ±0.5 stands for the remainder dependence on the number 
of flavors at 0{a^) for all physically reasonable A^ and is completely negligible for phe- 
nomenology. The last term is also very small ~ 0.4(as/7r)^. The running coupling can be 
parametrized in the standard form ( p.l7| ) with A~ = 1.41A^. 

Using the FAC approach (Grunberg, 1980, 1982, 1984), we rewrite eq. ( pl6|) as follows. 

..eff 



R{s) 



f 



1 ± ^ ± Oiat 

TT 



(9.7) 
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where the 3-loop correction is absorbed into the definition of the effective couphng given 
by eq. ( |2.17| ) with the A replaced by 



1 r2 



1.02A 



MS" 



As one can see, the new scheme MS almost coincides with the effective one and the fastest 
convergence is guaranteed within the wide range of energy defined by the renormalization 
group invariant criteria 

s s 
~ ^ 1 



MS 



The similar analyses can be done for the semi-hadronic decay j;ates of the r lepton 
calculated to 0{al) in section 7. The result for the ratio Rr in the MS scheme reads 



i?^ = 3(0.998 ±0.002) 



TT 



TT 



+ 9.83 



TT 



+ 0(at 



and to be compared to eq. (|7.10|) . Note that the as{Mz) is parametrized with the A 



(9. 



1"41Ams- 



MS 



In Fig. 14 we plot one-, two- and three-loop approximants to the T^j^f^^ in terms 
of the running quark mass (eqs. (|3.12|) - (|3.16| ), with = 5 and frif = nih) vs. the scale 
parameter t (Surguladze, 1994b). 
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FIG. 14. The approximants of the F^^.^ vs the scale parameter t 



One can see that the higher order corrections diminish the scale dependence from 40% 
to nearly 5%. The solid curve, corresponding to the three-loop result, became fiat in the 
wide range of the logarithmic scale parameter t. Moreover, the choice t = (MS-scheme) 
satisfies Stevenson's Principle of Minimal Sensitivity (Stevenson, 1981). 

Let us now try to estimate the theoretical uncertainty in calculations of R by the 
last included term in the_corresponding perturbative expansion. We get for the QCD 
contribution within the MS scheme the following result. 

«B - = ? + (v)" - ('^-^ ± (t)' * (^«- = 4%). (9.9) 

The analysis of Fig. 13 shows that the deviation of the four-loop approximant from the 
constant is also about 4% within a reasonably wide range of the t-parameter. This is 
consistent with Stevenson's principle. One should note that the above error estimate is 
only for the massless quark limit. There arc several difi^erent types of additional contri- 
butions, including those due to nonvanishing quark masses. This may change the above 
error estimate. All of the necessary information on the status of the additional corrections 
can be found in Kniehl (1994b, 1995b). 

As we have already mentioned, recently Brodsky and Lu (1994, 1995) have found the 
relations between the effective couplings cka and for the physical observables A and B 
in the following form. 

aA(AiA) = aB(AiB)(^l + Ta/b^ H ^- (9.10) 

The ratio of the scales of the corresponding processes /Ia/ /is is chosen according to the 
BLM scale setting prescription so that Ta/b is independent of the number of flavors. Thus, 
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evolving and Ob, they pass the quark thresholds at the same scale. It is shown that 
the relative scales satisfy the transitivity rule 

fJ'B fJ'C fJ'B 

So, C may correspond to any intermediate theoretical scheme such as MS, MS, etc. and 
the perturbative results can be tested without a reference to them. One of the impressive 
results of this method is a surprisingly simple relation between the effective couplings for 
the quantities R and Rr to the next-to-next leading order (Brodsky and Lu, 1994, 1995) 



TT 



IT 
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For more details and the relations between various other observables we refer to the original 
works by Brodsky and Lu (1994, 1995). 



On scheme invariant analyses 

Let us now outline the original method of scheme-invariant analyses for the perturbation 
theory results by Stevenson (1981a,b, 1982, 1984). We note first, that our analyses of 
perturbation series for R{s) and Rr has been done in the previous subsection within the one 
parametric family of the MS type schemes, where all (3 function coefficients are the same for 
any modification of MS. In the PMS method, renormalizarion scale and scheme dependence 
is parametrized by the scale parameter /i/A and the scheme dependent coefficients of the 
/3 function /?2; Ps, ■ ■ ■• Then the Principle of Minimal Sensitivity is applied to the variation 
of the above parameters and to 0{al) the "optimized" scheme corresponds to a flat two 
dimensional surface. Our curve for i?3 in Fig. 13 is just a one- dimensional slice at the 
particular MS value of the f32- The main points of the PMS formalism is as follows. (For 
the scheme invariant analyses of R{s) to 0{al) see Mattingly and Stevenson, 1994). To 
use familiar standard notation, we rewrite eq. (|2.7| ) for the couplant a = as(//)/7r 

dci 

6— = -6a^(l + ca + C2a^ H ), (9.11) 

where 

r = Mn^, 6 = 2/5o, c=^ (9.12) 

and for any modiflcation of the minimal subtraction prescription, the scheme dependent 
coefficient C2 = (32/(3o- The scheme and scale can now be parametrized by the quantities 
RS = (r, C2, C3, ...). The Principle of Minimal Sensitivity can be written as 

, = 0- (9.13) 
d(r;c2,C3, ...) 
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The number of scheme-scale parameters in the above equation is strongly correlated with 
n. Indeed, it is not difficult to show that the following self-consistency condition should 
hold for the nth approximant 

dRn ^/n+l\ 



0(a" 



(9.14) 



d{RS) 

This shows that the perturbative coefficients rj can depend on renormalization scheme 
only through parameters r; C2, Ci_i. Applying the Principle of Minimal Sensitivity in a 
form ( |9.13 ) to the approximants R2 and and taking into account ( 9.14 ), one finds that 
the quantities 



Pi 



r2, 

c 



(9.15) 



are renormalization scheme independent. Similar invariants can be constructed at each 
order of perturbation theory. The choice of r as a function of the ratio /i/A emphasizes 
that the renormalization scheme dependence involves only the ratio of these quantities 
and the optimization deals with r but not 11. The "optimal" values of renormalization 
scheme parameters r and C2 are defined by the following equations. To 0{a1), 



dRoir) 



To 0{al 



dr 



dRsi^T, C2] 



0. 



dr 
dRsir, C2) 



9C2 



(9.16) 



(9.17) 



(9.18) 



C2=C2 



Solving the above equations along with eqs. (|9.15|) for the renormalization scheme invari- 
ants and eq. (|9.11 ) for the couplant with the truncated MS (3 function, using the MS 
values of r2 and r^, one finds the "optimized" values of r, C2 and corresponding "opti- 
mized" approximants to O(a^). The theoretical error can be estimated, as in the previous 
subsection, by the last calculated term. One obtains the following "optimized" result for 
the QCD contribution in i?(34 GeV) in the massless quark limit (Mattingly and Stevenson, 



1994; Stevenson, 1994). 



rPMS 

"qcd 



0.051 ±0.001. 



(9.19) 



It is important to note that the above optimization procedure yields a negative value 
for the P2 invariant. This results in the existence of a solution of equation 



7 _ 

- + ca* + 3p2{a 



■*\2 







(9.20) 



with respect to a* - the value of the couplant for which the optimized third order (3 function 
vanishes. This allows, in principle, to do some analyses for R{s) at the low energies ^/s — > 
(Mattingly and Stevenson, 1992). 
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Finally, we also mention that the FAC approach (Grunberg, 1980, 1982, 1984) is a 
special case of the PMS (Stevenson, 1981a,b, 1982, 1984) method. Indeed, in the FAC 
approach all higher order approximants are equal to the effective couplant (compare to 
eqs. (p.l6|) and (|9.18|) ). From eqs. (|9.15|) one gets pi = r and p2 = C2 in the FAC approach. 



Conclusions 

In the present article we reviewed the current development of calculational meth- 
ods, algorithms and computer programs which allow one to evaluate the characteristics 
of the phenomenologically important physical processes to higher orders of perturbative 
QCD. We have considered Z hadrons, t~ ^ Ur + hadrons, H — > hadrons. The de- 
scribed methods are applicable to a wide class of calculational problems of modern high 
energy physics. We outlined the analytical three- and four-loop calculations for the above 
mentioned processes. 

The methods of analytical perturbative calculations available at present allow, in 
principle, one to evaluate various decay rates, cross-sections, coefficient functions in the 
operator product expansion, renormalization group functions etc. up to and including five- 
loop level. This would correspond, for instance, the decay rate in the process Z hadrons 
to O(a^). It seems that such a high order will completely fit the experimental state 
of the problem in the observable future. Indeed, for example, the 4% estimate of the 
theoretical error for the decay rate of Z-boson is based on the O(a^) calculation. The 
present experimental error at LEP is about 5%. 

The involvement of the heavier quarks in the physical processes makes it necessary to 
develop methods for calculation of the Feynman graphs with the propagators of massive 
particles. The expansion in terms of large or small masses may not always give satisfactory 
results. 

The problem of the renormalization group ambiguity of the perturbation theory re- 
sults and various methods for resummation of higher order corrections is a subject of 
growing interest and discussions in the literature. 

The future development of analytical programming tools towards the full automation 
of high order calculations would be welcome. This would greatly reduce the chance of errors 
in the calculations. On the other hand, the computer package with full implementation 
of the algorithm of high order analytical perturbative calculations would make it realistic 
to step up by one more order. 

We recognize that it is unavoidable that some of the relevant references have not been 
mentioned. We assure the reader that this is only due to our unintentional ignorance. 
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